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The Numerical Simulation of General Relativistic Shock Waves by a 
Locally Inertial Godunov Method Featuring Dynamical Time Dilation 

Abstract 

We introduce what we call a locally inertial Godunov method with dynamical time 
dilation, and use it to simulate a new one parameter family of general relativistic shock 
wave solutions of the Einstein equations for a perfect fluid. The forward time solutions 
resolve the secondary reflected wave (an incoming shock wave) in the SmoUcr- Temple shock 
wave model for an explosion into a static singular isothermal sphere. The backward time 
solutions indicate black hole formation from a smooth underlying solution via collapse 
associated with an incoming rarefaction wave. As far as we know, this is the first numerical 
simulation of a fluid dynamical shock wave in general relativity. 



-iv- 



Acknowledgments 

I would like to thank all the great teachers I have learned from throughout my academic 
career. They made many sacrifices for their students and rarely get the credit and respect 
they deserve. Without them, I would have fallen off this path of enlightenment and never 
written these words. A special thanks goes to my advisor Blake Temple. His insight, 
guidance, and enthusiasm was a real inspiration, and I feel very privileged to have worked 
with such a great man. I am also grateful to all my committee members for taking the time 
to look over my work. I want to send a very special thanks to the late Evelyn Silvia, whose 
passion for teaching was a big influence on my teaching endeavors. 

I want to send a warm thanks to the math department staff, especially Celia Davis and 
Perry Gee, whose efforts on making the department seem like a family did not go unnoticed. 

I wish to give huge thank you to my family, whose impact on my life over the years I 
am just beginning to understand and appreciate. I very special thanks goes to my mom, 
Kady Vogler. Her death was the most significant moment of my life, and it is a shame she 
missed one of my greatest accomplishments. May her soul rest in peace. 

Finally, I would like to thank to all my friends that have shown me love and support 
throughout the years. Their presence keep me sane throughout the challenging world of 
mathematics. Though there are many graduate students to thank, a special one goes out 
to Spiros Michalakis and Ernest Woei; many great ideas and conversations came out of 
those weekly lunches. My greatest appreciation goes to the Ricss household and the circle 
of friends that evolved from there. They are a second family to me, and they keep me 
grounded with simple pleasures. I want to pay homage to the game of basketball for giving 
me an avenue to exercise my muscles besides the brain. I want to thank all the bailers out 
there that influenced my game; this list includes, but not limited to, Duane Kouba, the 
math basketball team, and the ARC legends. 



-V- 



CHAPTER 1 

Introduction 

This thesis presents a numerical study of a canonical family of spherically symmetric 
solutions to the Einstein equations of General Relativity, such that shock wave formation 
is simulated in forward time and black hole formation is indicated in backward time. To 
start the simulations, we introduce a new one parameter family of initial data satisfying the 
constraints of the Einstein equations, and to perform the simulations, we develop a locally 
inertial Godunov method incorporating a new dynamical time dilation feature to account 
for the relativistic effects of curvature. 

Our point of departure is the work of Groah and Temple |GT04i| . who gave the first 
existence theory for shock wave solutions of the Einstein equations, starting from general 
initial data of bounded variation in the density and velocity. The setting is restricted to 
spherically symmetric spacetimes in standard Schwarzschild coordinates. For the analysis, 
Groah and Temple introduced the locally inertial Glimm scheme, a numerical method based 
on approximating spacetime by flat Minkowski space in each grid cell, such that the method 
was amenable to detailed mathematical estimates sufficient to prove convergence to weak 
(shock wave) solutions of the Einstein equations. The curvature of spacetime in the locally 
inertial Glimm scheme is accounted for by coordinate transformations between the locally 
inertial approximations in each grid cell. This technique is pedagogically interesting be- 
cause it parallels the so called correspondence principle, the physical principle that general 
relativity should reduce to special relativity in sufficiently small neighborhoods of space- 
time. Groah and Temple concludes in |GT04j that the curvature of spacetime necessarily 
becomes discontinuous at shock waves, resulting in solutions solving the Einstein equations 
in the weak sense of the theory of distributions. In this thesis, we develop these locally iner- 
tial formulations of the Einstein equations into a viable numerical method, what we refer to 
as a locally inertial Godunov method with dynamical time dilation, and we introduce a new 



one parameter family of interactive solutions on which we demonstrate the convergence of 
the method. 

Our numerical demonstrations are backed up by a general theorem which we prove. 
Namely, we prove if a sequence of approximate solutions, generated by our locally inertial 
Godunov method converges pointwise almost everywhere to a function along with the total 
variation of the fluid variables remain uniformly bounded under the limit, (exactly what 
we can demonstrate numerically), then the limit solution is an exact weak solution to the 
Einstein equations. By this general theorem, we need only demonstrate numerical conver- 
gence to a limit, with bounded oscillations, in order to conclude our simulated solutions 
accurately represent exact (weak) solutions of the Einstein equations. 

The one parameter family of spacetimes on which we test the convergence of the method 
are interesting in their own right. This one parameter family is constructed by matching 
initial data from a critical (k = 0) Friedmann-Robertson- Walker (FRW) spacetime to initial 
data for a Tolmann-Oppengeimer-Volkoff (TOV) spacetime at a point where the spacetime 
is only Lipschitz continuous, the smoothness required for our convergence theorem. Since 
our numerical method is based on standard Schwarzschild coordinates, this requires a coor- 
dinate transformation from FRW coordinates to standard Schwarzschild coordinates. Such 
a transformation requires an integrating factor, and we use the integrating factor derived by 
Smoller and Temple in |ST95] . (we call the resulting standard Schwarzschild coordinates 
version of FRW the FRW-1 spacetime). We also derive a new integrating factor that leads 
to a different standard Schwarzschild coordinates version of FRW, (which we call FRW-2). 
The initial discontinuity that separates FRW from TOV in the standard Schwarzschild co- 
ordinates initial data then generates a region of interaction between the two spacetimes, a 
region for which there is no closed form solution, and can only be constructed by numerical 
simulation. The ultimate goal of this thesis is then to simulate this region of interaction, 
and demonstrate convergence of the locally inertial Godunov method. 

The numerical simulation of the region of interaction between FRW and TOV is a prob- 
lem motivated by the original work of Smoller and Temple in |ST95j , further developed in 
|ST03|, IST04) . In [ST95j , the first exact shock wave solution of the Einstein equations is 
constructed by matching an FRW spacetime to a TOV spacetime across a shock interface, 
such that the resulting matched solution, given by exact formulas, is a true weak solution 



of the Einstein equations. The resulting matched spacetime (we denote as the FRW/TOV 
metric) with a shock wave was a simple model for an explosion into a static singular isother- 
mal sphere, relevant to models of star formation jST95j . Based on Smoller and Temple's 
work, we built a computer visualization [JSV] to demonstrate the qualitative behavior of 
their model. The Smoller and Temple (Sm/Te) construction requires taking a different 
sound speed on the FRW and TOV sides of the shock wave finely tuned to reduce the 
region of interaction between FRW and TOV to a single pure shock wave. As a result, they 
only required the existence of, not an explicit formula for, the integrating factor. In our 
simulation, we take the same equation of state on both sides, a more realistic assumption, 
and as a consequence, a region of interaction, for which there are no exact formulas, between 
the lightlike curves emanating from the initial point of connection between FRW and TOV 
metrics is generated. Based on this background, we interpret our numerical FRW/TOV 
model simulation in forward time as resolving the secondary reflected wave in the Sm/Te 
model of an explosion into a static singular isothermal sphere. In particular, our simulation 
confirms for the first time the refiected wave is a shock wave, not a rarefaction wave. We also 
simulated the time reversal of the FRW/TOV solution. In this case our simulation indicates 
the formation of a black hole in the resulting solution, which consists entirely of rarefaction 
waves. Based on our simulation, one can only expect the observed strong rarefaction wave 
evolving toward the coordinate center will form a black hole. Note that a black hole cannot 
be computed exactly in standard Schwarzschild coordinates because black holes generate 
coordinate singularities in this coordinate system. A proposal for future work would be to 
simulate the black hole all the way by transforming the solution to Eddington-Finkelstien 
or Kruskal coordinates [Wei72) . 

We start our discussion by briefly introducing the Einstein equations in Chapter [2l In 
this introduction, we formulate these equations in the simplest setting for shock waves in 
General Relativity; this setting is the spherically symmetric Einstein equations with the 
stress energy tensor for a perfect fluid and an equation of state p = ap, p = const. In this 
setting, the Einstein equations reduce to a conservation law with source terms coupled with 
a system of ODEs. 

Chapter [3]sets out the background information for the FRW, the TOV, and the FRW/TOV 
matched models. In particular, we construct a one parameter family of initial data that 



agrees with the FRW and TOV metrics in standard Schwarzschild coordinates on either side 
of a single point where they match continuously. The continuity of the metric at the match- 
ing point implies the initial data meets the constraints of the Einstein equations introduced 
by Groah and Temple in [JGB06] . More specifically, the initial data is built by connecting 
a FRW metric with a TOV metric across an initial discontinuity in the fluid variables, such 
that the metric components are matched Lipschitz continuously. This one parameter family 
of initial data enables us to construct the one parameter family of interacting spacetimes 
that solve the Einstein equations. To be simulated by our locally inertial Godunov scheme, 
the metric must be in standard Schwarzschild coordinates, and while the TOV is already 
in this coordinate system, the FRW metric needs to be mapped over. We demonstrate 
two mappings of the FRW metric over to standard Schwarzschild coordinates, denoted as 
FRW-1 and FRW-2. With these two coordinate transformations, the one parameter fam- 
ily of shock wave solutions, the FRW/TOV metric, has two different forms based on the 
coordinate system used. These two forms are used as a further test of the correctness of 
our numerical simulations, by verifying the simulations produce the same spacetime in two 
different coordinate systems. 

In Chapter [H we introduce the locally inertial Godunov method. This chapter begins 
with a discussion of the solution to the Riemann problem |Smo83| for relativistic compress- 
ible Euler equations in Minkowski spacetime. The Riemann problem is the fundamental 
building block of the Godunov method, |LeV92j . We outline the solution to this Riemann 
problem as demonstrated by Smoller and Temple in |ST93j . In order to incorporate the 
Riemann problem in the Godunov method for simulating non-flat spacetimes, we augment 
the solution to the Riemann problem with the time dilation feature. To this end, we de- 
termine the effects of time dilation on the Godunov step and incorporate it into our locally 
inertial Godunov method. The chapter ends with a layout of the locally inertial Godunov 
method and all the steps necessary to implement it. This method is a fractional Godunov 
scheme that defines the solution inductively and contains four major steps: a Riemann 
problem step, a Godunov step (with time dilation), an ODE step, and an update step. 

In Chapter [5l we prove the main theorem that establishes the consistency of our locally 
inertial Godunov method with time dilation. It states, assuming the convergence of our 
approximate solution constructed using our method with a total variation bound at each 



time step, the convergent solution is indeed a weak solution to the Einstein equations. The 
two main difficulties of this proof are handling the jump in the fluid variables across each 
time step and showing the jump in the metric across each space step is accounted for by the 
term added to the ODE step. The assumptions to this theorem are obtained numerically by 
the simulation results. Thus, our theorem is perfectly tailored to our numerical simulation: 
we numerically establish convergence and a total variation bound, and this main theorem 
proves that assuming these conditions, we can conclude convergence to a weak solution of 
the Einstein equations. 

In Chapter [6l we begin the numerical study of the locally inertial Godunov method with 
dynamic time dilation by simulating the pure non-interacting FRW and TOV spacetimes 
individually. Each of these continuous models provide us with a different scenario to test 
and fine tune our method. These models are also important for correctly handling the 
non-interaction regions and the ghost cells for the subsequent shock wave model. Within 
this chapter are graphs of the solutions and tables of the convergence results, along with 
the explanations behind them. 

In Chapter [71 the first glimpse of shock wave solutions to the Einstein equations is 
shown by simulating the FRW-l/TOV model. We demonstrate this one parameter family 
converges and produces quantitatively different solutions as we vary the free parameter. 
Regardless of the parameter chosen, the solution always contains two shock waves, an 
incoming and outgoing wave, enclosing a region of higher density, denoted as the interaction 
region. Over and above the convergence and graphical results, we show the interaction 
region converges to the cone of sound (the extent of the sound like information emanating 
from the initial discontinuity), and we also show outside this interaction region the FRW- 
1 and TOV metrics are preserved under the simulation. Changing coordinate systems, 
the FRW-2/T0V model simulation is similar to the FRW-l/TOV model. We discover and 
numerically confirm the FRW-2/T0V model is the same solution as the FRW-1 /TOV model 
that differs by a non-linear coordinate transformation, and this second model provides a 
pedagogically interesting numerical confirmation of the covariance of the Einstein equations 
in standard Schwarzschild coordinates. 

In Chapter [HI we examine the time reversed FRW/TOV model. The numerical simula- 
tion shows this solution consists of two rarefaction waves, one outgoing and one incoming. 



encompassing a region of lower density. In particular, the simulation indicates there are no 
shock waves present. Therefore, this simulation shows the underlying solution is smooth 
and thus a strong solution to the Einstein equations, unlike the forward model. We record 
the convergence of this solution and the convergence of the interaction region to the cone of 
sound along with the preservation of the non-interaction regions. The simulation convinces 
us the incoming rarefaction wave will form a black hole in finite proper time. That is, 
we extend the time of the simulation up until extreme relativistic effects of time dilation, 
characteristic of a black hole, appear. We cannot simulate all the way up to black formation 
in standard Schwarzschild coordinates because black holes are coordinate singularities in 
standard Schwarzschild coordinates |ST97b| , but the metric component A is monotonically 
decreasing to a small value before time dilation makes our numerical calculations imprac- 
tical. Indeed, we numerically demonstrate the mass function near the black hole formation 
gets within 1.084 < 1.125 = 9/8 of the Schwarzschild radius of that mass; thus, we have 
numerically confirmed the solution gets within the Buchdal statiblity limit of 9/8ths the 
Schwarzschild radius. Buchdal's theorem states that whenever the mass of the star gets 
within 9/8ths of the Schwarzschild radius for that mass, there is no static configuration 
with a finite pressure at the center capable of holding the mass up |ST97bj . The impli- 
cation then is a star that gets within 9/8ths the Schwarzschild radius cannot be prevented 
from collapsing into a black hole; thus, we view our results here as a numerical physical 
proof that black hole formation must occur in the time reversed FRW/TOV model. Thus, 
the simulation indicates a black hole forming out of a smooth solution to the Einstein equa- 
tions. With the success of running the forward and reverse FRW/TOV models, we have 
simulated general relativistic analogs of both the 1 and 2 family of shock and rarefaction 
waves, that is, the analogs of all the elementary waves of conservation laws. This fact gives 
us confidence in the ability and accuracy of our method to simulate any solution to the 
Einstein equations for a perfect fiuid in standard Schwarzschild coordinates. 

In Chapter [9l we put dimensions back into the problem, giving our simulations of the 
earlier chapters a physical context. For convenance and concreteness, we set the Einstein's 
speed of light constant and Newton's gravitational constant both to one. To place units on 
our numerical values requires us to understand the effect this has on our units. The effect of 
this is the unification of the units of time, length, and mass into one set of units, which we 



choose as the units of mass. Understanding this unification allows us to recover the other 
units, time and length. To end this chapter, we consider our simulation on the solar and 
galactic scale by transforming our numerical values to familiar units, giving us a firm grasp 
of the scale and magnitude involved within our simulations. 

In Appendix[Al we outline the code used to perform the simulations in this paper. A CD 
is attached that contains all this code, which is approximately 8,000 lines. The appendix is 
designed to provide the reader an overview of the organization of the code and to help the 
reader navigate to specific areas of interest. 
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CHAPTER 2 

Preliminaries 

This chapter introduces the background to the topic of this thesis, the locally inertial 
Godunov method with dynamical time dilation. This is a fractional step numerical method 
that approximates the convective part of the fluid dynamics in each grid cell by the solution 
of a Riemann problem in a (flat) Minkowski background spacetime, such that each cell is 
endowed with a time dilation factor. Averages of the fluid variables are taken, and a second 
fractional step then evolves these averages according to evolution by a dynamical system 
that accounts for the undifferentiated source terms. The main issue is at shock waves the 
gravitational metric has a jump discontinuity in the derivative, so the curvature of the 
gravitational metric becomes discontinuous at shocks, and the purpose of this thesis is to 
confirm our method is effective at numerically simulating general relativistic shock waves. 

The locally inertial Godunov method with dynamical time dilation is based on the locally 
inertial formulation of the Einstein equations for spherically symmetric spacetime metrics in 
standard Schwarzschild coordinates and was first derived by Groah and Temple in [GTOOj . 
Groah and Temple provided an exposition of their work along with the work of Smoller 
and Temple in |JGB06] . All the equations in this chapter are taken from the joint work 
of Groah, Smoller, and Temple in |JGB06] . In this study, we consider the stress energy 
tensor a perfect fluid with an equation of state p = a^, where a = Const, the simplest 
possible setting for shock wave propagation in general relativity. In this case, the particle 
number density and entropy decouple from the energy and momentum equations, and the 
fluid equations close at the level of the Einstein equations alone. This equation of state 
models flow at constant temperature, and in the case a = it models the important 

cases of free particles in the extreme relativistic limit, as well as the case of pure radiation, 
c.f. [Wei72| . In |JGB06] . Groah and Temple show that the four standard Schwarzschild 
coordinates Einstein equations are weakly equivalent to the system of PDE's obtained by 

^C.f. |JGB06] where is used in place of a 
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taking the two equations that express the vanishing divergence of the stresses in place of 
the two equations that involve time derivatives of the metric. This results in a system of 
conservation laws with source terms in which the special relativistic Minkowski variables 
are taken as the conserved quantities. Fortunately, the equations close because all time 
derivatives of metric cancel out for this particular case. The beauty of these equations is 
that they reflect the locally inertial character of spacetime. The locally inertial Godunov 
method with dynamical time dilation is then a numerical method tailored to the simulation of 
solutions of these equations. We prove the consistency of this method, and demonstrate its 
value as an efficient method for computing shock wave solutions of the Einstein equations. 

In Section [2.11 we introduce the Einstein equations, and we proceed to show how these 
equations with spherical symmetry and a simple equation of state reduces to a system of 
PDEs which express the locally inertial character in Einstein's theory. In Section 12.21 we 
formulate the conservation law with source terms. Finally, we pose the general problem 
that is conducive to shock wave solutions to the Einstein equations, the starting point of 
our studies. 



The fundamental equations of general relativity, the Einstein field equations, is given 



where G is the Einstein curvature tensor, T is stress energy tensor (the source of gravity), 
and K is Einstein's coupling constant 



where Q is Newton's gravitational constant, and c the speed of light. For concreteness and 
convenience, we take the speed of light as c = 1 and Newton's gravitational constant as 
Q = 1 throughout this paper. The consequences of setting these constants to one is discussed 
in detail in Chapter [9j These equations describe how the mass/energy of an object curves 
space and how the curvature, in turn, stretches or squeezes matter in the three spatial 
directions. 



2.1. Introduction to Einstein Equations 



by 
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The component form of the Einstein equations is 
(2.3) Gij{x) = KTij{x), 

where 

(2-4) Gij = Rf^j - -RZlQij 

denotes the individual components of the Einstein curvature tensor. We let i, j G {0, 1, 2, 3} 
refer to the individual components in a given coordinate system. Throughout, the Einstein 
summation convention is used where repeated up-down indices are summed from to 3, 
unless otherwise noted. 

Modeling the sources of the stress energy tensor by a perfect fluid, T is given by 

(2.5) T'^ = {p + p)w'w^ + pg'\ 

where w denotes the unit 4-velocity vector of the fluid (tangent vector to the world line 
of the fluid particle), p denotes the energy density (measured in the inertial frame moving 
along with the fluid), and p denotes the fluid pressure. 

For simplification, the only metrics under consideration are spacetime metrics g which 
are spherically symmetric. A general spherically symmetric metric takes the form 

(2.6) ds"^ = gijdx'dx^ = -A{t, r)dt'^ + B{t, r)dr'^ + 2D{t, r)dtdr + C(t, r)dn'^ , 

where the metric components A, B, C, D are functions of the time and space variables (t, r) 
only, dQ"^ = dO"^ + sin^ 6d(j)'^ represents the standard line element of the unit 2-sphere, and 
X = {x^, . . . ,x^) = {t,r,6,(j)) is the spacetime coordinate system. A spherically symmetric 
metric is in standard Schwarzschild coordinateqj if it is written in the more simple form of 



(2.7) ds^ = -B(t, r)dt^ + -r^-^dr"^ + r^dn"^. 



A classic argument transforming a spherically symmetric metric (j2.6p into standard Schwarzschild 
coordinates (12. 7p , which is thoroughly explained in Chapter [3l 



^Note that the choice of notation in |JGB06] for standard Schwazschild coordinates was ds^ — —A{t, r)dt'^ + 
B{t,r)dr^ +r^dn'^ 
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In this paper, the simplest setting for shock wave propagation in General Relativity is 
assumed: a spherically symmetric metric in standard Schwarzschild coordinates ()2.7p where 
the sources are modeled by a perfect fluid (|2.5|) with an equation of state 

(2.8) P = crp, < fj < c, 

where a is assumed to be constant and y/a represents the sound speed. Using MAPLE to 
put the standard metric (12. 7p into the Einstein equations (j2.1j) gives us the following system 
of four coupled partial differential equations, 

(2.9) ^{-4 + l-l} = .i.^T00, 



(2.12) ^ i rA + B" -<t \ =2^r AT^^ 

where the quantity <I> is 

2i AB 1 (aV B' BA' B fB'V BB'A' 

(2.13) '^-A^^^A^~^[a) "T"7T+2"V^y "2"5"T' 

where "prime" represents ^ and "dot" represents ^. 

The mass function is M = M(t, r) is defined through the identity 

2M'^ 



(2.14) A= [1 



r 



which is interpreted as the total mass inside radius r at time t. Using this total mass, 
equation ()2.14p leads to an equivalent form for equation (j2.9p given by 

(2.15) M' = Inr^BT^^, 



2.1. INTRODUCTION TO EINSTEIN EQUATIONS 



12 



and equation (j2.10p becomes 

(2.16) M = -^Kr^BT^\ 

Assuming a perfect fluid (j2.5p . the components of the stress energy tensor T*-' satisfy 
the following relations 

(2.17) = Ito? 



(2.18) T'' = \l^n} 



(2.19) T^^ = ATI} 

(2.20) = ^, 

where represent the components of the stress energy tensor in flat Minkowski spacetime. 
When an equation of state takes the form ()2.8p . the components of Tjv/ are given as 

('2 211 ^00 _ c + q- ^ 



(2.22) = 



— 



(2.23) rl) = ^^cV. 



We denote the fluid velocity v, in place of the 4- velocity vector w, as measured by 
an observer fixed with respect to the radial coordinate r. Combining (|2.15p together with 
()2.17p . a formula for the mass is 

/"T 

(2.24) M{t, r) = M{t, tq) + ^ / T^^{t, rydr. 

(Equations (|2.9p - (|2.24p are taken from equations (1.3.2)-(1.3.15) in the work of Groah, 
Smoller, and Temple |JGB06] .) We note this equation ()2.24p combined with ()2.14p is an 
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easier way to compute the metric component A than using equation ()2.9p . and these equa- 
tions are used to compute A in our locally inertial Godunov scheme presented in Chapter 

m 



2.2. Statement of the General Problem 

The presence of shock waves in solutions of (j2.9|) - (|2.12|) makes the stress-energy tensor T 
discontinuous, and consequently, the metric components A and B are Lipschitz continuous 
at best. Since ()2.12p contains second derivatives of A and B, this equation can only be 
satisfied in the weak sense. In |JGB06] . Groah and Temple show when the metric is 
Lipschitz and the stress-energy tensor is bounded in sup- norm, system (j2.9p - (j2.12p is weakly 
equivalent to a system in which equations (I2.10p and ()2.12p are replaced by 

2 



-01 
M ) 



(2.25) {T'^n,o + {^T'J} ^ = 

{rS},o+{VlBTl,i}^^ = 

Here, {}^o find represent the derivative with respect to the time and space variable, 
respectively. The other two equations ()2.9p and ()2.1ip are rearranged as 



(2.27) ^ = - ^T-, 



B' (i - 1) Kx n 

(2.28) + 

The variable x is written instead of r when the equations are expressed as a system of 
conservation laws consistent with the literature on the subject. This conservation form is an 
advantageous formulation because equations (j2.25p and ()2.26p form a system of conservation 
laws with source terms, where the conserved quantities u = (T^j,T^j) are independent of 
the metric (c.f. (j2.2ip and (j2.22p ). Moreover, all the time derivatives of the metric (j2.10p 
and (|2.12p cancels out, and the source terms to the conservation laws (j2.25p and ()2.26p 
are independent of any derivatives, resulting in a solvable system. The resulting system 
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(j2.25p - (|2.28p can be written compactly as 

ut+f{A,u)x = g{A,u,x), 

(2.29) 

A' = h{h.,u,x). 

Here, the conserved quantities are the Minkowski energy and momentum densities 

(2.30) n = (TOO,rOi) = (nO,^z^), 
and the metric components are written as 

(2.31) A = {A,B). 
The flux function becomes 



(2.32) f{A,u) = VABiTlj,T}^ 



M )■ 

The source terms of the conservation laws, (j2.25p and ()2.26p . are 

(2.33) g{A, u, x) = (5° (A, u, x),g' (A, u,x)), 
where 

(2.34) g'>iA,u,x) = --VABTij, 

x 

and 

g^{A,u,x) = 

(2-35) ^ 2kx 

2 I X X A 

The LHS of the metric ODEs, (12:271) and <!(I?M . are 

(2.36) h{A,u,x) = {h (A, u, x), h^{A,u,x)), 
with 

(2.37) /i°(A, u, x) = - ^txT??, 

X 
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and 



(2.38) 



h {A,u,x) = — 



B ({l-A) 



X 



+ KxTi 



ill 

M 



} 



Thus, for weak solutions to the Einstein equations, the coupled PDEs (I2.9p - (l2.12jl are equiv- 
alent to a conservation law with source terms paired with a system of ODEs (j2.33p . In this 
thesis the convergence of the locally inertial Godunov method (defined in Chapter H]) is 
shown to be a consistent and effective first order method for computing solutions of system 
()2.33p . equivalent to the Einstein equations when shock waves are present and the gravi- 
tational metric is C^'^. Such solutions only solve the Einstein equations in the weak sense 
of the theory of distributions. Our general system (|2.33p applies to the case of spherically 
symmetric metrics in standard Schwarzschild coordinates, and in the next chapter, a list of 
candidates meeting this criteria is assembled. 
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CHAPTER 3 

One Parameter Family of Shock Wave Solutions 

This chapter presents the theory behind the one parameter family of shock wave so- 
lutions of the Einstein equations whose numerical simulation is the subject of this thesis. 
Our point of departure of this chapter is the paper [ ST95| in which Smoller and Temple 
(Sm/Te) construct exact spherically symmetric shock wave solutions to the Einstein equa- 
tions for a perfect fluid. Our purpose here is to include a more realistic equation of state, 
and thereby simulate the secondary reflected wave, which requires numerical simulation be- 
cause it cannot be expressed in closed form. The Sm/Te solutions were realized by matching 
two distinct metrics to form a Lipschitz continuous hybrid metric |ST97a|, IJGT03] . In 
particular, the critical (k=0) Friedmann- Robertson- Walker (FRW) metric for the expanding 
universe was matched to the Tolman-Oppenheimer-Volkoff (TOV) metric, a model of the 
interior of a star, across a shock wave interface. To perform this matching, one is required 
to provide a coordinate map of a spherically symmetric metric, like the FRW metric, over 
to standard Schwarzschild coordinates. Recall, a general spherically symmetric metric takes 
the form 

(3.1) ds^ = -A{t, r)dt^ + B{t, r)dr^ + 2D{t, r)dtdr + C(t, r)dJ^^ 

where the metric components A, B, C, D are functions of the time and space variables (t, r) 
only. The expression dil^ = d9^ + sin^ Odcj)'^ represents the standard line element of the unit 
2-sphere. A metric in standard Schwarzschild coordinates takes the form 

(3.2) ds^ = -B(t, f)dP + df^ + f^dn^, 

where again the metric components are functions of the time and space alone. Throughout 
this paper, the barred coordinates {t, f) are reserved for metrics in a standard Schwarzschild 
coordinates system (13. 2p while the unbarred coordinates (t, r) are for metrics that are only 
spherically symmetric (j3.ip . This FRW metric transformation in not only necessary for the 
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matching, but it is also a prerequisite for the locally inertial Godunov scheme since it is only 
capable of simulating metrics in standard Schwarzschild coordinates. When constructing 
these shock wave solutions, Sm/Te proved the existence of the necessary coordinate trans- 
formation, but they did not require any detailed information about it. Here, we do not have 
this luxury since this information is essential to run the simulation. Our main goal of this 
chapter is to uncover these details on transforming the spherically symmetric FRW metric 
into standard Schwarzschild coordinates. We discover there exists at least two of these 
transformations of the FRW metric, denoted as FRW-1 and FRW-2. In the subsequent 
chapters, these FRW metrics expressed in different standard Schwarzschild coordinates are 
used to test our locally inertial Godnuov method, by simulating these continuous metrics. 

In Section [3?H we briefly discuss the FRW metric. The FRW-1 metric is obtained by pro- 
viding a candidate for the transformation of the FRW metric into standard Schwarzschild 
coordinates and proving the validity of this transformation. The fluid variables also get 
mapped over as functions of the new coordinates. In Section 13. 2^ we review the standard 
procedure for mapping a spherically symmetric metric to standard Schwarzschild coordi- 
nates, and the FRW-2 metric is constructed using this procedure. The other metric in the 
matching, the TOV metric, is discussed in Section [3.31 This metric is already in standard 
Schwarzschild coordinates, but the main aspects of it needed in this paper are explained in 
this section. We conclude with the details of matching the FRW and TOV metrics contin- 
uously across a shock wave interface in Section 13. 4[ It is in this section, we make precise 
the one parameter family of shock wave solutions by defining the shock wave surface. 



3.1. The FRW-1 Metric 

The conformally flat {k = 0) FRW metric takes the form 

(3.3) ds^ = -dt^ + R\t){dr^ + r^dn'^}. 

The fluid density, as a function of time t, associated with this metric takes the form 
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where again 

(3.5) K = 

is Einstein's coupling constant (j2.2p . To determine the fluid velocity, the fluid in this 
metric is comoving relative to the background. More precisely, its velocity vector, w = 
{w^ ,w^), has the form = 0, for i = 1, 2, 3, or, in words, the movement of the fluid 
is limited to the time coordinate and none of the space coordinates. The FRW metric g 
being diagonal and the velocity vector having unit length implies 

(3.6) «;0 = ^/3^=l. 

This velocity vector is expressed in terms of the classical coordinate radial velocity v = dr/dt 
of the particle paths of the fluid as f = 0. 

Recently, Smoller and Temple derived a new set of equations with solutions describing a 
two parameter family of expanding wave solutions of the Einstein equations which includes 
the critical Friedmann universe in the standard model of cosmology [STJ. In deriving this 
set of equations, they discover a coordinate map for the FRW metric (j3.3p over to standard 
Schwarzschild coordinates (j3.2p for cr = |, and where R{t) = \/i is the scale factor. We 
state their result and derivation in the following theorem 



Theorem 3.1.1. Assume we have an equation of state p = and k = 0. Then the 
FRW metric /i3. 3]) under the coordinate transformation 



r 



--Vtr, 



(3.7) _ . _2 ^ ^2 



4t^ j 4 

goes over to the following metric in standard Schwarzschild coordinates 

(3.8) ds^ = ^df + — ^df2 + r'^dn'^, 

1 — f 1 — 

where the fluid velocity v is related to f/t by 
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Proof. To map the FRW metric tensor (j3.3p from the unbarred coordinates {t,r) over 
to the barred coordinates (t , f ) , we use the tensor transformation law 

(3-10) ~'-^ = '^^dE^dii^- 

More specifically, at each point, g transforms by the matrix transformation law 
(3.11) g = J'^gJ 

for a bilinear form since the Jacobian matrix J = transforms the vector components of 
the x-basis {gf^} into their corresponding coordinates of the x-basis {^}- The inverse 
of the Jacobian matrix is easily computed by taking partial derivatives of the coordinate 
transformation equations (|3.7p . resulting in 



(3.12) J 



So the Jacobian matrix is 
where 

(3.14) \J-^ 



At - 



The metric in barred coordinates g is computed as 



Vi -TT^ \ / -1 \ / Vi 



r 



(3.15) 

1 I -it-i) 

|J-1|2 1 

Using (|3.14p to solve for 
(3.16) \J- 



L 1 / \ t / \ ^- 1 



t-^ 



.l,2_(i-T)^ 
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gives us 



(3.17) 



V 







T 

^ 4t 



Note, the new metric is a function of r/t, and it is useful to define a new variable rj as 
(3.18) 



_ r \/tr r 
As a function of rj, the metric becomes 



(3.19) 







V 







1-^ / 

or equivalently it is written in standard Schwarzschild coordinates form (13. 2p . 



(3.20) 

with 

(3.21) 



-^dP + — ^-df2 + f'^dn'^ 



A{t,r) = 1 



B{i,r) 



^ — nl' 

^ 4 



The metric is also a function of the fluid velocity v, shown by finding the relationship 
between v and our variable rj. To this end, the old velocity vector w is related to the new 
velocity vector w using the tensor transformation law 



(3.22) 



w 



dx' 



-w 



so for our comoving velocity vector w = (1,0,0,0) this law simplifies to 

dx" 



(3.23) 



w 



Using ()3.23p with our coordinate transformation (13. 7p . the new velocity vector coordinates 
are 



(3.24) 



w 







dx^ dt 
dx^ dt 



= 1, 



_i dx^ dr 1 r _2 _q 

w = — -TT = — = w = w = 0. 

dx^ dt 2^t 
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For a spherically symmetric metric in standard Schwarzschild coordinates (j3.2p the radial 
velocity v = dr/dt is computed as 

Hj^ 1 1 r If 

(3.25) I 



VAB 2 2 t • 

For our metric p.20p . A = , implying V AB = 1. By (|3.18p . the relationship between v 
and rj is 

(3.26) ^ = |. 

We substitute this relation (j3.26p into ()3.20p to prove ()3.8p . By defining the variable 
^ = f/i, the coordinate transformation (j3.7p implies 

(3.27) 



t (l + g)t 1 + ^2' 
proving (|3.9p . □ 



The FRW metric in standard Schwarzschild coordinates (13.2 



(3.28) ds^ = + — ^df2 + f^dO^ 

1 — t;^ 1 — 

with 

1 



(3.29) A{0 = l-v{i)\ B{0 



is a function of ^ alone. We consider now the fluid variables {p,v). The fluid velocity v is 
implicitly a function of ^ by equation (13. 9[) . but through manipulation the fluid velocity u 
is a function of ^ explicitly 



(3.30) viO = ^ 

where the minus sign is taken to ensure v < 1. The density p is a function of the FRW time 
t (j3.4p and is transformed into a function of by multiplying the density function (j3.4p by 
to obtain 

(3 31) , 3^_3^_MO! 
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using (j3.26p . As a side note, this equation along with (j2.15p imphes that the function ^ 
at a point, where M is the total mass, can also be written as a function of ^. Rearranging 
equation (|3.31|) . the density function becomes 

Mi? 



(3.32) p{iJ) 



These results are recorded as a Corollary to Theorem 13.1.11 

Corollary 3.1.1. The fluid variables {p,v) corresponding to the FRW metric in stan- 
dard Schwarzschild coordinates \3. ^) are given as 

(3-33) p{S.,r) = — =2~' ^(0 = h ' 

where ^ = f/t. 

3.2. The FRW-2 Metric 

In this section, we show the existence of another mapping of the FRW metric (13. 3p 
into standard Schwarzschild coordinates. As opposed to the a priori knowledge of the 
FRW-1 transformation in Theorem I3.1.H the FRW-2 transformation is developed through 
construction. There is a classic argument in |Wei72| to build a coordinate transforma- 
tion (t, r) — )• (t, f) taking a general spherically symmetric metric (13. ip over to standard 
Schwarzschild coordinates p.2p . and this argument is repeated to build the FRW-2 trans- 
formation. In particular, we take the spherically symmetric FRW metric 

(3.34) ds'^ = -dt^ + R^{t){dr^ + r^dQ^}, 
over to the standard Schwarzschild coordinate form 

(3.35) ds'^ = -Bit, f)dP + , dr^ + r'^dn'^. 
^ ^ ^ A(t,r) 

We first match the spheres of symmetry 

(3.36) f2fifi2 ^ R^r^d^l^, 

so the first component of the coordinate mapping (t, r) — )• (t, f) is defined by the following 



(3.37) 



r = r[t, r) = R(t)r. 
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Rewriting ()3.37p as 

(3.38) r = ^, 

we differentiate and square to obtain the following differential form relationships, 

1 Rr 

(3.39) dr = —df - —rdt 

R R^ 

and 

(3.40) dr^ = -^df^ - '^drdt + ^^dt'^ 
Substituting (I3.40p into ()3.34p . results in the FRW metric in {t,f) coordinates 

(3.41) = - (^1 - dt^ + dr^ - ^dtdr + f'dn\ 

In order to complete the coordinate transformation (t, r) {t,f), t = t{t, r) must be defined 
to eliminate the cross term dtdf in (j3.4ip . Since it is already matched, the f^dJl term is 
ignored in our subsequent calculations to simplify notation. We show the procedure for 
eliminating the cross term for a generalized metric of the form 

(3.42) ds'^ = -C{t, f)dt^ + D{t, f)dr'^ + 2E{t, f)dtdr, 

and handle our specific case afterwards. To eliminate the cross term in p.42p . an integrating 
factor ^ = ^{t,f) is chosen to satisfy the equation 

(3.43) |(^^) = -|(*^)' 
implying 

(3.44) dt = "^{Cdt - Edf] 

is an exact differential, proven at the end of this section in Lemma 13.2.21 Exactness 
gives us the existence of t{t,f), completing the definition of our coordinate transformation 
(t,r) — )• (t, f). We consider options for the integrating factor ^{t,f) and the corresponding 
coordinate transform t = t{t, r) later; for now, we continue mapping the metric over, by 
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assuming satisfies ()3.43p . Squaring both sides of (|3.44p results in 

(3.45) dP = ^^{C^dt^ - 2CEdtdr + E^df^}, 
and rearranging leads us to 

1 E'^ 

(3.46) - Cdt^ + 2Edtdr = --^^r^df + ^df^. 

Substituting this equation into ()3.42p gives us the diagonal metric in the barred coordinate 
frame 

(3.47) d~s' = -^df + + ^) dr'. 
For our specific case, we have 

i?^f^ Rf 

(3.48) c=l-^, D = l, E = - — , 

and substituting these into the general form (j3.47p gives us 

(3.49) ds'^ = ^df + ^df^ + r^dn^. 

^2 i_R^] I i?^r-^ 



For the pure radiation phase, where a = 1/3, the cosmological scale function becomes 
R{t) = \/t, and the resulting metric is 

(3.50) ds^ = ^ ^df + — i — ^df2 + f^dVL^. 

The definition of r] (j3.18p allows us to rewrite the metric in the barred coordinates as 

(3.51) ds^ = -^d? + , , df^ + f^d^^. 



Notice this metric (j3.5ip closely resembles the FRW-1 case (j3.20p except for inclusion of 
^, indicating our ability to recover the FRW-1 metric by choosing the correct integrating 
factor. With this in mind, we look for solutions of the integrating factor equation (|3.43p . 
In light of the functions C{t,f) and E{t,f) defined in (|3.48p . equation (j3.43p becomes 
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The following lemma expresses the existence of two solutions for this integrating factor 
equation 



Lemma 3.2.1. The following PDE for the function '^{t, r) 

has at least two solutions, a constant solution ^'i(t,f) = ^'q and a dynamical solution 
(3.54) ^2(t,f) = 1'o. 



4*2 + f2 ' 

where is some constant. 



Proof. The constant function, r) = ^'o, solves the PDE (|3.53p since 



To find the other solution for the integrating factor ^'2(^,7'), a simplification of (|3.53|) is 
obtained by noticing 



d f \ d 



r 



so ^'2 satisfying (13.53P is equivalent being a solution to the following PDE 

^'■'^^ -dFV-^^)=^KYt 

We proceed by constructing a function ^'2 satisfying ()3.57p . Suppose ^'2 has the form 

(3.58) ^2(t,r-) = ^, 

Vr 

where f{rj) is an unknown function of the predefined variable rj = f/t. Taking the partial 
derivatives of 13. 581 leads to 

^ ' ' df tVf 2fi' dt ^ \ t^J ' 
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where the prime represents differentiation with respect to the variable ry. Substituting rj 
into equation (j3.57p gives us 

Plugging in the partials (j3.59p into the PDE (|3.60p and simplifying results in the following 
ODE 

f'iri) 4-772 

(3.61) J \ u _ I 



fill) 277(4 + r?2) 
Using partial fractions, the solution of this ODE is 



(3.62) /(r/) = /01 



'4 + 7/2 ^ 



r/o 4 + 7/2 

Combining the constants together by defining 



(3.63) *o = /oi 

produces the unknown function /(?/) as 



,2 



(3.64) /(77) = ^o^^. 

Therefore, the second integrating factor becomes 



(3.65) ^I/2(t,r-) = M = xi7 / /? vl;, / * 



f(4 + 7/2) V 4*2 + f2 ' 

proving (|3.54p . □ 
The first integrating factor ^'i, the constant solution, gives us a differential form (j3.44p 

of 

T'2 \ 7" 



(3.66) dt = ^!{Cdt - Edf} = | (^1 - dt - —df 

Since this is an exact 1-form, the coordinate transformation t(t, f) is obtained through 
integration as 

(3.67) t{t,f) = -^Jl + ^\t, 
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or 

(3.68) i(t,r) = ^'o(l + Qi, 

where f = R{t)r = \/tr is used to connect the two. By setting = !> we recover the 
FRW-1 coordinate transformation (j3.7p into standard Schwarzschild coordinates, providing 
another proof for its existence. 

The second integrating factor ^'2, the dynamical solution, has the corresponding 1-form 
(pmi) of 

(3.69) £ = *{« - E<lf) = *„ { (1 - 5) * - (^) *} ■ 

Again, since this 1-form is exact, there exists a function t{t, f) such that 



, „ , dt ^ / t f dt ^ / t 

(3.70) =,1, /^_^ 1- T- = -*c 



5t V + f2 V 4*2 y 5f V 4*2 + f2 \2t 

To solve for t we integrate the ^ equation to obtain 



(3.71) t = ^^^^+g{t), 

for some function g{t). Taking ^ of (j3.7ip leads to g{t) = C, the constant function. By 
setting g{t) = 0, the coordinate transform corresponding to the second integrating factor 
(1331)1 is 



(3.72) f(,.) = :^^if±i!. 

Hence, our new coordinate transformation of the FRW metric (|3.3|) into standard Schwarzschild 
coordinates is 

f =Vtr, 

(3.73) 



- ^-0 /4t2 + f2 

Similar to the FRW-1 case, there exists a relationship between v and r/. The velocity 
vector transformation laws, (j3.22p and (j3.23p . along with our new coordinate transformation 
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(j3.73p provides the new velocity vector w coordinates 



(3.74) = — = — = ■ — 

The radial velocity v = dr/dt is computed as 



dr 1 r 



w 



w 



(3.75) 



w 



1 



V4t + 



1 



V 



4*2 + f2 



1 



1 



'AB 2^QVt ^/AB 2^0 V t y/AB 2* ^AB 2 ' 

where we used (|3.5ip that \/AB = 1/^*. Remarkably, this relationship matches the FRW-1 
case (I3.26p . Using this relationship, the FRW-2 metric (I3.5ip is rewritten in terms of the v 
as 



(3.76) 



ds' 



1 



-dP + 



1 



1 



rdf2 + f^dQ^ 



^2(1 _ ^2) 

Unlike the FRW-1 metric, the fluid variables are not a function of the ratio ^ = f/t, so 
the FRW-2 metric relies on the unbarred coordinate time t. The coordinate transformation 
p.72p is used to find t as a function of the new coordinates (t , f ) , 



(3.77) 



t{i,f) 



t2 + y|rr72^ 



All these results for the FRW-2 metric are recorded in the following thereom 



Theorem 3.2.1. Assume we have an equation of state p = and k = 0. Then the 
FRW metric \3. 3|) under the coordinate transformation 

f =\/tr, 

2 V i 

associated with the integrating factor 



(3.79) ^'(t,f) = *o ' * 



4t2 + f2 ' 

goes over to the following metric in standard Schwarzschild coordinates 

(3.80) ds^ = -- r^TT ^dt^ + T^dr^ + ^^dn\ 

^ ' ^'2(1- 1)2) l--y2 
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The fluid variables {p, v) corresponding to this metric are 

3 



(3.81) 



where the unbarred time coordinate is the following function of {t, r) 
(3.82) 



2*J 



This section ends with the lemma stating any ^' satisfying (|3.43p makes dt (I3.44p into 
an exact 1-form. 



Lemma 3.2.2. If the integrating factor ^{t,r) satisfies 
(3.83) 

then the differential from dt = {^C)dt — {^E)df is exact. 



|(*C) + |(*E)=0, 



Proof. Suppose ^ satisfies (I3.83p . By the following computation 



(3.84) 



curliy) 



i j k 

dt df dz 



|(*c:)-|(-*E)|k = o 



this assumption is equivalent to ^ satisfying curl{v) = 0, where v = (^C, —"^EjO) is a 
vector in 3-space. To simplify notation, define functions g{t,f) = ^E and h{t,f) = —^E. 
Since curl{v) = 0, there exists a function f{t,f) such that 



(3.85) 



df df 

-^ = ^C = g and = -^E = h. 

at or 



Let C be a curve in {t,f) space parameterized by a, denoted by C = {t{a),f{a)), going 
from the point Pi = (t(ai), f(ai)) to the point P2 = {t{a2),r{a2)), shown in Figure [3Tl 
Now take the line integral of the 1-form dt along the curve C to obtain 



(3.86) 



rP2 ra2 
dt= gdt + hdf = 

C J Pi Jai 



g{t{a),f{a))^ + h{t{a),f{a))^ 



da 



Using the chain rule along with (I3.85p . we have 

(3.87) dt = £ ^da = f{t{ai),r{ai))) - /(^(aa), r(a2))). 
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t 




r 



Figure 3.1. Arbitrary curve in (t, r)-space 



Hence, di is only dependent on the endpoints of the curve C which means di is inde- 
pendent of path; therefore, the 1-form dt is closed. Since in subsequent chapters we only 
consider a convex domain D = {rmin ^ x < rmax^t > 0}, Poincare's Lemma [Rud76] 



3.3. The TOV Metric 

For the outer TOV solution in our shock wave simulation, we use the general relativistic 
static isothermal sphere with the equation of state p = for the pure radiation phase 
as derived by Smoller and Temple in |ST95j . For completeness we summerize the results 
required here. This TOV metric has the form 



states a closed 1-form on a convex domain is exact, proving the claim. 



□ 



(3.88) 




The time metric component B takes the form 



(3.89) 



B{r) = Bo{r)^ 



where ^/a is the speed of sound, and the mass M function is given as 



(3.90) 



M(r) = 4TT^r, 



where the parameter 7 is a constant dependent on the equation of state constant. 
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which agrees with equation (3.4) in |ST95| . The fluid variables in the TOV metric are 

(3.92) P(^) = ^' ^ = 0' 

where the velocity is zero since the TOV metric is static. 

Unlike the FRW metric ()3.3p . the TOV metric is already in standard Schwarzschild 
coordinates (13.21) with 



(3.93) A{r) = l-^^^^, B{r) = Bo{f)^. 



both independent of the time coordinate t. Notice since 
(3.94) 2£M(r) ^ ^^^^^ 



this equality simplifies the A metric component to be 
(3.95) A{r) = 1 - 871^7, 

a constant, independent of f, for constant a. 



3.4. One Parameter Family of Shock Wave Solutions 

With both the FRW and TOV metrics in standard Schwarzschild coordinates, we dis- 
cuss the matching that stitches the two together. For notation purposes, we use Apjiw and 
BpRw to denote metric components of the FRW metric in standard Schwarzschild coordi- 
nates for the general form ()3.5ip . regardless of which integrating factor is chosen, and Atov 
and Btov to denote the corresponding components for the TOV metric (j3.88p . Assume the 
(t, f) coordinates, describing the TOV and FRW metrics in standard Schwarzschild coordi- 
nates, represent a single coordinate system where both metrics are matched continuously 
across a shock surface (i.e. ApRw = Atov and BpRw = Btov)- Matching the A metric 
component gives us 

(3.96) Afrw = l-v^ = l — = Atov- 
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Recall for both FRW metrics, FRW-1 and FRW-2, the relationship v = r/2t holds. Using 
this relationship along with the density equation (j3.4p . the square of the velocity becomes 



(3.97) = ^pf^ 



Substituting (j3.97p into (j3.96p and rearranging produces 



(3.98) M = ^pr\ 



By taking c = 1 and using equation ()3.5p to replace Av, this equation simplifies to 

(3.99) M(f) = —p{t)f^. 

3 

The equation (j3.99p defines the shock surface f = f(t) implicitly in the (t, f) coordinates. 
To obtain the corresponding curve in the original FRW coordinates {t,r)^ the substitution 
of f = R{t)r into (j3.99p is made to define the shock surface r = r[t). Interestingly, equation 
(j3.99p is independent of the integrating factor meaning the shock surface is the same 
curve in (t, f) space regardless of the integrating factor chosen. 

To obtain the other metric component must be matched at the shock surface f[t), 

(3.100) Bfrw = .T,2,/ 2^ = ^o(r)^ = Btov- 

This matching allows us to finish the definition of ^ by defining the integrating factor 
constant ^'o once the TOV time scale factor Bq is chosen. In the simulation setup, we 
reverse the roles of these constants; we choose ^'o such that the coordinate speed of light 
(i.e. \/AB) is one on the FRW side, for both FRW-1 and FRW-2, and equation (i3T00]l 
defines Bq. Notice by choosing \I'o, which defines the rate at which time progresses uniformly 
across the matched metric, there remains one free parameter f or t related by the shock 
surface equation (j3.99p . It is in this sense there exists a one parameter family of shock 
wave solutions. In this paper, we consider the initial position of the discontinuity r as our 
parameter. With the equations laid out for the FRW-1, FRW-2, TOV, and the matched 
metrics in standard Schwarzschild coordinates, we posses the required information for our 
locally inertial Godunov method, discussed in the next chapter. 
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CHAPTER 4 

Locally Inertial Godunov Method 

This chapter is dedicated to the algorithm for the locally inertial Godunov method 
featuring dynamical time dilation. Although in this thesis, the method is used to simulate 
the family of shock waves identified in Chapter [3l the method can be applied to simulate 
general spherically symmetric flows. This method is a modification of the locally inertial 
Glimm method by Groah and Temple |JGB06l, IGT04) , and many of the equations within 
this chapter are taken from their work. In their paper, Groah and Temple devise the 
fractional Glimm method to prove the existence of shock wave solutions to the spherically 
symmetric Einstein equations for a perfect fiuid. This method is a technique to evolve 
solutions from initial profiles for the conserved quantities u{t, x) and the metric A(t, x) 
in standard Schwarzschild coordinates satisfying the Einstein equations. Recall, standard 
Schwarzschild coordinates take the form 



(4.1) ds^ = -B(t, f)dP + , dr^ + r^dn"^. 

A{t,f) 

Here, we modify the locally inertial Glimm method by replacing the random choice Glimm 
step in Groah and Temple's work with an averaging Godunov step to obtain more consistent 
and less jagged solutions, with the ultimate goal of simulating these shock wave solutions 
in General Relativity. Also, the purpose of the Groah and Temple construction is to prove 
an existence theorem; therefore, the development lacks some details needed to numerically 
construct the solution. A goal of this chapter is thus to fill in these details, like an algorithm 
to find the middle state to the Riemann problem. Our method can be interpreted as a locally 
inertial scheme, in the sense it exploits the locally fiat character of spacetime. Each grid 
cell is considered a locally fiat frame, and we handle the time dilation between frames, by 
choosing a reference frame relative to which time can be synchronized. The frame chosen 
is the one in which the factor ^/ AB is one. We denote this reference frame as a unitary 
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frame. We also denote the factor V AB as the coordinate speed of light to represent the rate 
at which hght travels relative to the unitary frame. Throughout this paper, we will refer 
to it as the speed of light, not to be confused with the speed of light constant c = 1, and 
it will be clear from context which one is being used. Note that for any single frame, the 
metric component B (j4.ip can be rescaled by a change in the time coordinate to make the 
(coordinate) speed of light one in that frame, so the reference frame is quite arbitrary and 
not so important in tracking time. However, it is important what this reference frame's 
speed of light is relative to the other frames around it, causing time dilation between the 
frames. These ideas will be explored and expanded throughout the chapter. 

In Section [4. II we discuss the initial value problem in Special Relativity. In this section, 
we study solving the Riemann problem in Minkowski spacetime, which is a unitary frame. 
Solving the Riemann problem is at the heart of the Godunov step of our method, and 
its importance cannot be overstated. Section 14.21 explains time dilation and how it affects 
the averaging of our Godunov step, extending our ability to perform this step on non- 
unitary frames. With the background material set, we conclude in Section 14.31 by stating 
the algorithm of the locally inertial Godunov method. This method formulates the solution 
inductively and has four major steps: a Riemann problem step, a Godunov step (with time 
dilation), an ODE step, and an update step. 

4.1. The Initial Value Problem in Special Relativity 

In this section we develop solutions to the relativistic compressible Euler equations in 
flat Minkowski spacetime for the case p = ap with a constant a 




together with the initial conditions 

(4.3) piO,x) = po{x), v{0,x) = vo{x). 



Note, this equation (j4.2p is our conservation law ( ()2.25p and (j2.26p ) without the source term 
where y/AB = 1. This problem is a specific case of the initial value problem for a general 
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system of nonlinear hyperbolic conservation laws in the sense of Lax [Smo83j , 



ut + {F{u))^ = 0, 



(4.4) 



■u(0, x) = uq{x). 



In our case, we have 



(4.5) 




+ 1], p{a + c^)^^] 
c — V J 



and 



(4.6) 



F{u) ^ F') = (p{a + c^)^^, p[{a + c')^ 



To distinguish between the two sets of variables, we refer to the pair (p, v) as the fluid 
variables and to the pair {u^,u^) in (j4.5p as the conserved quantities. Both of these variables 
play an important role in the locally inertial Godunov method. The conserved quantities 
are needed in implementing the Godunov step, and the fluid variables are needed to solve 
the Riemann problem along with giving us more physical meaning behind our simulation. 
Fortunately, Groah and Temple showed [ JGB06] there is a 1 - 1 correspondence between 
the fluid variables and the conserved quantities, giving us the following result 

Proposition 4.1.1. The mapping {p,v) — t- {u^,u^) is 1 - 1, and the Jacobian determi- 
nant of this mapping is both continuous and non-zero in the region p > 0, \v\ < c. 

In the locally locally inertial Godunov method, we constantly transfer back and forth 
between the conserved quantities and the fluid variables, so the inversion of the mapping 
{p,v) — >• {u^,u^) defined in (|4.5p is necessary to find the fluid variables as a function of the 
conserved quantities, stated in the following corollary to Proposition 14.1.1] 

Corollary 4.1.1. The mapping {u^,u^) — t- {p,v) takes the form 



(4.7) 




{{a + c^)u^ 



,0 



V(cT + c2)2(uO)2 -4a(ni)2}. 



(4.8) 



p{u^u') 



(c^ — v'^)u^ 
(cr + c'^)v 
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Proof. Using the definition of the conserved quantities ()4.5p . both are solved as equa- 
tions for p 

\ I P / ,9\ 9 o-2i2' 

and 

(4.10) f = K^- 

Setting these equations, (|4.9|) and (|4.1U|) . equal to each other and simplifying gives us the 
following quadratic equation in v 



(4.11) ^u\'^ - (a + c^)u^v + c^u^ = 0, 



so there exists two candidates for v 

(4.12) V = -^{{a + c2)n0 ± J{a + c2)2(u0)2 - 4cj(ni)2}. 

In order to determine which solution to use, we normalize the speed of light (i.e. c = 1) 
and rewrite (I4.12P as 



(4.13) „.i!!(^±i) 4("')'^'^ 



2u^a y Y (u0)2(a + l)2y 

Since u'^ > and 0" + 1 > it, the following inequalities hold " 2u^f7^ ^ ^ ^^"^ ^ 
In order to keep v less than the speed of light (i.e. v < 1), the minus sign is taken in ()4.12p 
to obtain (j4.7p . Substituting this v into (|4.10p gives us p as (|4.8p . □ 



The main focus of this section is to solve the Riemann problem for system (|4.2p . The 
Riemann problem is the initial value problem with initial data uq{x) of a pair of constant 
states separated by a jump discontinuity at x = 0, 



(4.14) 



uo{x) 



UL X < 

un X > 0. 
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Note, in light of Proposition I4.1.1t the conserved quantities ul and ur are uniquely deter- 
mined by the fluid variables {pl^vl) and {pr,vr), and the symbol u can be interpreted as 
either set of variables throughout this paper. 

A general theorem of Lax |Smo83| states for any system of conservation laws (j4.4p 
which is strictly hyperbolic and genuinely nonlinear in each characteristic field, the Riemann 
problem for this system has a unique solution in the class of elementary waves if ur and ul 
are sufficiently close. The proof for this result is a constructive one, relying on the structure 
of the state space. More specifically, given a point ul in u-space, there exists a family of 
i-wave curves connecting the point ul to any other sufficiently close point ur by traversing 
a 1-wave curve followed by a 2-wave curve, where the middle state um is the intersection of 
these two curves. Depending on the direction along the 1-wave curve, a shock or rarefaction 
wave is between the states ul and um, and similarly for the 2-wave curve, the direction 
taken determines whether a shock or rarefaction wave is between the states um and ur. As 
a side remark, each i-wave curve has second order contact at the point ul (i-e. their first two 
derivatives are equal at this point) |Smo83j . Figure ITT] gives us a sample of the web-like 
structure for these family of curves for a 2 system of conservation laws referred to as the 
(non-relativistic) p-system |Smo83] . for p = 1. In this figure, the red dot represents the left 
state Ul. The 1-shock curve is represented by the blue/red graph, and the 1-rarefaction is 
represented by the red/cyan graph. The 2-shock and 2-rarefaction curves, originating from 
the 1-wave curve, are represented by green and yellow graphs, respectively, with the red 
stripped curves representing the ones emanating from the left state. One can interpret the 
coordinate system of wave curves defined for each ul as giving us a road map providing 
directions from our starting point ul to our destination ur: the 1-wave curve is the first 
street traveled and 2-wave curve next street traveled, and the middle state um is where 
the streets traveled intersect. The solution becomes the left state connected by a 1-wave 
to the middle state which is connected by a 2-wave to the right state. This result is only 
a local one because for certain conservation laws (including 7-law gases), not all the points 
in u-space can be connected by this web structure, due to the formation of the so-called 
"vacuum" states |Smo83j . 

For other systems, this local connectivity of points in u-space can be extended to a 
global connectivity of points, where the "vacuum" states do not appear. Fortunately, for 
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our system of interest (j4.2p . when p = ap, the "vacuum" states do not appear and our 
i-wave curves cover the entire u-space. Following Smoller and Temple |ST93j . Groah and 
Temple extend in |JGB06] this general theorem of Lax our system 



Theorem 4.1.1. There exists a solution of the Riemann problem for system with 
an equation of state p = ap, < ^/a < c, as long as ul and satisfy 

(4.15) PL > 0, pR> 0, 
and 

(4.16) — c < vl < c, — c < vr < c. 
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Moreover, the solution is given by a 1-wave followed by a 2-wave, satisfies p > 0, and all 
speeds are bounded by c. This solution is unique in the class of rarefaction and admissible 
shock waves. 



The main goal of this section is to use the details within the proof of this result to 
construct the explicit solution to the Riemann problem for any left and right states. This 
section is organized into three subsections. Subsection 14.1.11 covers preliminary material 
needed to solve the Riemann problem. Finding the middle state along with the type of 
elementary waves connecting it to the left and right states is discussed in Subsection 14.1.21 
Bringing all these details together, Subsection 14.1.31 explains how to solve the Riemann 
problem at any point of interest. 



4.1.1. Preliminaries. In this subsection, equations needed to solve our Riemann prob- 
lem ()4.2p and (14.140 are presented. Along with the fluid variables, there is another set of 
variables, called the Riemann invariants, important in solving the Riemann problem effi- 
ciently. The Riemann invariants are functions which are constant along the integral curves 
or rarefaction curves in the state space |Smo83j . The Riemann invariants r and s for our 
system (|4.2p are 



(4.17) r(p, v) = lln -^^ ln(p). 



(4.18) s{p, v)=^-ln (^^) + Hp), 
where 

(4.19) K 



2ac 



,2 



(CJ + C2)2- 

There is a typo in the statement of the Riemann invariants in the Groah and Temple paper 
|JGB06j in equations (2.5.73), (2.5.74), (4.2.12), and (4.2.13). As well as moving freely 
between the conserved quantities and the fluid variables, the ability to go back an forth 
between the fluid variables and the Riemann invariants is needed within the locally inertial 
Godunov method. 
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The fluid variables (p, v) are recovered from the Riemann invariants (r, s) by algebraic 
manipulation of (I4.17[) and (j4.18p . More specifically, subtracting (|4.17p from (|4.18p to solve 
for the density 



(4.20) p{r, s) = exp 



s 



'2K 

and adding (j4.17p to (j4.18p to solve for the velocity 

c(l - 6^'+'' 



(4.21) v{r,s) 



1 + 6^+'' 



To solve the Riemann problem, it is necessary to know the speed of the shock waves 
in the solution to determine which side of the discontinuity the point of interest is located. 
The speed for the 1-shock is the following function of beta (3 



(4.22) si 
while the speed for the 2-shock is 

(4.23) S2 



'/+(/g) + ^ 
/+(/?) + 



'/-(/3) + ;S 



V/-(/3) + ^ 

These shock speeds, (j4.22p and (j4.23p . are calculated in a frame where the particle velocity 
V is zero. To obtain these quantities in an arbitrary frame, the Lorentz transformation law 
for velocities must be applied. The transformation law is given as follows: if in a Lorentz 
transformation, the barred frame moves with velocity v measured in the unbarred frame 
with Si as the speed of the i-shock wave measured in the barred frame, and if s denotes the 
speed of the shock in the unbarred frame, then 

(4-24) . = 

The eigenvalues Xi{p,v) of the system (14. 2p are used to determine the speed of the 
rarefaction waves. These eigenvalues are Ai = —y/a and A2 = y/o' when the particle 
velocity is zero. Using the Lorentz transformation law (j4.24p for velocity v, the eigenvalues 
become 
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and 



(4.26) A2 



v + y/a 

I _|_ v^f 



To find the solution of the Riemann problem within a rarefaction wave, we need to solve 
for the fluid variables based on the eigenvalues and the Riemann invariants. More precisely, 
formulas for v as a, function of the eigenvalues, p as a function of (s, f ), and p as a function 
of (r, v) are necessary to solve for a state in a rarefaction wave. A quick calculation on 
(j4.25p and (|4.26p shows f is a function of the first eigenvalue 

(4.27) "(^0 = ^, 

or is a function of the second eigenvalue 

A2 - 



(4.28) v{\ 



2) 



1 



\f5\2 



Another quick calculation on (|4.17p and (|4.18p provides us with /o as a function of (r, v) 

{ pi ( 1 ( c + V 
(4.29) p{r, v) = exp { -\ -r; { r - - In 



and /9 as a function of (s, v) 



K \ 2 \c-v 



2 / 1 f c + f 
(4.30) p{s,v) = exp { J —{ s - -In 



K \ 2 [c-v 

4.1.2. Finding the Middle State. This subsection describes a numerical algorithm 
to find the middle state and the associated waves to the Riemann problem for our system 
(|4.2p . To accomplish this task, not only does one need to determine in which region the right 
state is located, governing the set of curve equations to use, but one also needs to solve these 
curve equations which cannot be solved explicitly. With this in mind, finding this middle 
state poses quite a challenge with the fluid variables {p,v). Fortunately, this problem is 
easier in the Riemann invariant coordinate system or the rs-plane. This coordinate system 
simplifles the i-wave curves in the state space because the rarefaction curves become straight 
lines, resulting in a clearer segregation of the state space. Another advantage is the i-wave 
curves are geometrically invariant across the rs-plane, so the shape of the wave curves are 
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independent of the base point {vltSl), implying an i-wave curve at a point in this plane can 
be mapped by rigid translation onto any other point. Within this subsection, we deal with 
the Riemann invariants as our variables, keeping in mind the goal is to obtain the middle 
state in the fluid variables. 

To distinguish between the different coordinate systems, we denote the fluid variables 
(/9, v) by lower case u and the Riemann invariants (r, s) by upper case U . More specifically, 
define 

Ur = {rR{UR),SR{UR)) = irR{pR,VR),SR{pR,VR)) 

(4.31) 

Ul = {rL{uL),SL{uL)) = {rL{PL,VL),SL{pL,VL)), 

where the transformations (j4.17p and (j4.18p are used. In the rs-plane, the 1-shock curve Si 
for the system is given by the following parametrization with respect to the /3, < /3 < oo: 

Ar = r-rL = -l M/+(2i^/5)} - Mln{Um = 5[(/3), 
(4.32) ^ 

A, = , - = -iln{/+(2K/3)} + y^ln{/+(/3)} = 5J(/3), 
and the 2-shock curve 5*2 is given by: 

Ar = r-rL = -\ ln{/+(2i^/3)} - ln{/_(/3)} ^ S'm, 

As = s-SL = -^ln{/+(2K/3)} + y|ln{/_(/3)} ee 5|(/3), 



(4.33) 



where 



(4.34) /^(/3)^l + /3|l=Fyi + | 
and 

(4.35) /3 = /3(.,..)^(^ + ^')' 



2a2 (c2_t,2)(c2_^;2)- 

In this coordinate system, both rarefaction curves are straight lines parallel to the coor- 
dinate axises along the positive directions in the rs-plane. More precisely, the 1-rarefaction 
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curve Ri is: 

Ar = r — ri = /3, 

(4.36) 

As = s — sl = 0, 
and the 2-rarefaction curve i?2 is given by: 

Ar = r — ri = 

(4.37) 

As = s - SL = /3, 

for the parameter /3, < /3 < oo. Notice within the equations for the shock and rarefac- 
tion curves (|4.32p - (j4.37p the differences Ar and As along the curves only depend on the 
parameter /3, proving the geometric invariance mentioned earher. For notation purposes, 
we define the following sets based at a point Ul for the shock curves 

(4.38) Si{UL) = {(r, s) : 3/3 > such that s = sl + 5f (/?), r = rL + S[ (/?)}, i = 1, 2, 

and the rarefaction curves 

Ri{Ul) ={(^ s) : s = SL and r > ri], 

(4.39) 

R2{Ul) ={(^ s) : r = rz, and s > sl]- 

Using this knowledge, a simpler family of curves is built to solve for the middle state [/m, 
as seen in Figure W?I\ 

Using these curves, the rs-plane is segregated into four disjoint open regions I, II, III, 
and IV as depicted in Figure HiS) In each region, a possible right state is placed and denoted 
by [/* for i = 1, 2, 3, 4. This figure also shows the 1-wave and 2- wave curves emanating from 
Ul as the thick lines with the 2-wave curves associated with the points C/* for i = 1,2,3,4 
emanating from the 1-wave curve as dashed lines. Region I corresponds to connecting the 
point Ul to the point [/^ by 1-rarefaction wave followed by a 2-shock wave. Region II 
corresponds to the l-shock/2-shock wave case, while region III and region IV corresponds 
to l-shock/2-rarefaction and l-rarefaction/2-rarefaction wave cases, respectively. Since the 
wave curves are geometrically invariant, we only consider the change from the right to left 
states, or the quantities Ar = r/j — r/, and As = sr — sl where the right state is U^ for 
i = 1,2,3,4, depending on the region under consideration. The middle state Um for all 
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Figure 4.2. Family of curves for our system in the rs-plane 



these cases is found by considering each quadrant in the rs-plane separately. After the 
quadrant breakdown, we discuss the technique used to find the correct /3(s) needed to solve 
the shock wave equations (j4.32|) and (j4.33p . Note that the shock curves cannot be solved 
exactly, and a threshold of e must be predefined, with the goal of solving the shock curve 
equations with a maximum error of e. 

We start with the point If^ in +/+ quadrant, which is determined by Ar > and 
As > 0. This quadrant is the same set as region IV, the l-rarefaction/2-rarefaction case. 
Since the rarefaction curves are straight lines, the middle state becomes rM = and 

SM = SL- 

Next we turn our attention to the point in the -/+ quadrant, where Ar < and 
As > 0. This quadrant is a subset of the region III, the l-shock/2-rarefaction case. Our 
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r 



Figure 4.3. Labeling the rs-plane 

goal is to find /? such that {rR,s*) £ Si{Ul)- After finding /3, the middle state will be 
irM,SM) = irR,s*). 

The point in the +/- quadrant, where Ar > and As < 0, is similar to in 
the -/+ quadrant. This quadrant is a subset of the region I, the l-rarefaction/2-shock 
case. Due to the geometric invariance of the 2-wave curve, we find a f3 such that the point 
{f*,sji) G S2{Ul), and the middle state will be {rM, sm) = {fL + {fji — r*), sl). 

Finally, consider the point f7^ in the -/- quadrant, where Ar < and As < 0. This 
quadrant is a superset of region II, the l-shock/2-shock case, so a mechanism is needed to 
determine when C/^ is not in region II but in region I or III as discussed below in the two 
shock algorithm. If C/^ is in region III, we find such that Um = {sm,''']\j) G Si{Ul) and 

such that (sj?,r/j) G S2{Um)- The middle state Um is, the intersection of the 1-shock 
curve Si{Ul) and the 2-shock curve S2{Um) crossing the right state Ur. 

When solving one shock curve equation, as in the +/- or -/+ quadrant, a standard 
bisection method |BF97] is used to find /3. We show how we implement the bisection method 
for the point f/^, where a 1-shock is followed by a 2-rarefaction. In this case, the bisection 
method is used to find the /3 that satisfies the equation Ar(/3) = — {ri + S'[(/3)) = 0, 
and our goal is to find a (3 such that Ar(/3) < e. The bisection method requires a starting 
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interval [/3mim /^max] where there exists a /3 within this interval such that Ar(/?) = 0. The 
shock curves, S'-{f3) and Sf{f3) for i = 1,2, are monotone decreasing functions of the variable 
/3. To find our interval, an initial guess, like /3 = 10^, is chosen and Ar(/3) is computed. If 
Ar(/3) > 0, our guess is too small (i.e. rji > ri + 5[ (/?)), and /3 is decreased. The power of 
our initial guess (3 = 10^ is decreased (i.e. k = 4,3, 2, 1, 0, —1, . . .) until a /c is found such 
that Ar(lO^) < 0, which means the guess is too big now. With this k, we set Pmax = lO'^'^"^ 
and /3min = 0. On the other hand, if for the initial guess Ar(/3) < 0, it is too big (i.e. 
fR<fL + 'S'i(/3)), and /3 needs to be increased. The power of our initial guess /3 = lO'^ is 
increased (i.e. k = 6, 7, 8, 9, . . .) until k is found such that Ar(lO^) > 0, which means our 
guess is too small now. With this k, we set [^max = and /3mjn = 10^~^. Either way, 
with the interval [f3max , l^min] established, the bisection method is implemented until a /3 
is found where Ar(/3) < e. For the -/+ quadrant, the same algorithm is used, but for the 
equation As(/3) = sr- {sl + S^{/3)) = 0. 

For the point U'^ in the -/- quadrant, the pair (/3^,/3^) solving the equations 

(4.40) Ar(/3S/32) - (r^ + S[{p') + S',{f)) = 
and 

(4.41) As{p\f3^) ^SR- (sL + StW') + 5|(/32)) = 

is sought out, and our goal is to find such that Ar(/3^,/3^) < e and As(/3^,/3^) < 

e. To perform the bisection algorithm for both (3^ and /3^, intervals [/5mm'/^max] ^^'^ 
[l^mini f^max] ^'^^ needed, just like solving for one shock equation. Since the equation (j4.40p is 
dominated by the parameter in order to get our interval [P^in, l^max]^ '^^ assume 
simplifying (|4.40p into Ar(/3^) = tr — (r^ + S\[f3^)) = and repeat the same procedure as 
above, increasing or decreasing the power k in our guess lO'^, for the parameter /3^. Simi- 
larly, we repeat this process to find [/3^m'/^max]' by setting Z?-*^ = in (j4.4ip and working 
with the equation As(/3^) = sr~ (s/, + 5'|(/3'^)) = 0. If either initial guess lO'^ gets too small, 
like k < —20, then we assume the point C/^ is in the wrong region. Depending on which 
parameter it is, /3i or (32, the U'^ must be in region I or region III, respectively. In either 
case, our problem boils down to only solving one shock curve equation, and it is handled 
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according to the above procedure. Equipped with our initial guess (/3^,/3^), we compare 
the errors Ar(/3\/32) and As(/3\/32) against each other. If Ar{/3^,l3^) < As(/3\/32), we 
perform a bisection method step on /3^; otherwise, we perform the step on /3^. We repeat 
this process until Ar(/3\/3^) < e and As(/3\/3^) < e. 

4.1.3. Solving the Riemann Problem. Being able to find the middle state with the 
connecting waves enables us to build the solution to our Riemann problem in its entirety. 
In particular, for an arbitrary point {t,x), the solution u{t,x) can be determined for any 
Riemann problem. This subsection covers the details of finding this solution at an arbitrary 
point. We also use this process to build a Riemann problem simulator. 

Suppose we have a Riemann problem (|4.14p to our conservation law (j4.2p and a point 
{t,x), and we want to solve for the state u* = u{t,x) = (p*,^;*) as a solution to this 
problem. Converting ul and ur over to Riemann invariants, using ()4.17p and (j4.18p . and 
implementing the above procedure, the middle state Um is determined along with the 1- 
wave connecting to the left state and the 2-wave connecting to the right state. Converting 
Um from Riemann invariant variables to fluid variables using ()4.20p and (I4.2ip gives us um- 

Without loss of generality, we consider the l-shock/2-rarefaction case, where the other 
cases are handled similarly. Since the solution to the Riemann problem is self-similar, 
the ratio of our point of interest s = x/t is compared to the wave speeds, usually done 
left to right. We first compare this ratio against the 1-shock speed si, computed by the 
transformation law ()4.24p with v = vl and the speed (sj) computed by ()4.22p using the /3 
found in the middle state algorithm. If s < si, then the state proceeds the 1-shock, and 
"U* = liL is the solution. Otherwise, we compute the speeds for the left and right sides of 
the 2-rarefaction fan. The left speed S2 is computed by (j4.26p with um as the input, and 
the right speed is computed similarly using ur instead. If si < s < , then the state 
is between the shock and rarefaction waves, and u^, = um- If < s, then the state is 
beyond the rarefaction wave, and = ur. li S2 < s < s^, then the state is within the 
rarefaction wave, and more work needs to be done to find the solution. We use the fact that 
within the rarefaction wave the solution varies smoothly, and every state u between um 
and ur moves at the speed \2{u) |Smo83| . Since the speed for our state is A2(ti*) = x/t, 
the inversion of the second eigenvalue (j4.28p gives us the fluid velocity v^,. With the fluid 
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Figure 4.4. The Riemann Problem Solver 



velocity, the density p^, is found by (I4.29P where r/j is used since this Riemann invariant 
is unchanged across the 2-rarefaction wave. For the rarefaction wave case, our solution 
becomes = (p*,^*). This completes the procedure to find the solution to the Riemann 
problem for any arbitrary point {t,x) for the l-shock/2-rarefaction case. 

With all this knowledge of the relativistic compressible Euler equations, we construct a 
Riemann problem simulator to test our solver and show us solutions to the Riemann problem 
for this system of equations in Minkowski spacetime. A glimpse of this simulator is shown in 
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Figure 1331 for concreteness, the 1-shock and 2-rarefaction case is displayed, where the left 
and right states are chosen to be ul = {pl,vl) = (10^,0.3) and ur = {pr,vr) = (10^,0.6), 
respectively. All the information needed to solve this Riemann problem is embedded into 
this figure. The colored rectangle is the spacetime cell containing the Riemann problem. 
The X-axis is horizontal while the t-axis is vertical, and the white half circle at the bottom 
center represents the origin {t,x) = (0,0). The speed of light is chosen to be one (i.e. 
c = 1) so twice as much space relative to time is needed to enable light to travel in both 
directions. The left/right state is displayed under the space time cell on the left/right side, 
and the middle state um = (2.002 x 10^,0.0639) is recorded above this cell. The density is 
displayed using a relative color map, where the highest value for the density (10^) is red, 
and the lowest value (10*) is yellow, with other values linearly interpolated between these 
two colors. The velocity is represented by the brown arrows where the length of the stem 
represents magnitude, and the direction left/right represents negative/positive velocity. The 
1-shock wave is represented by the blue line with the corresponding speed (si = —0.1717) 
recorded on the top left. The 2-rarefaction wave is represented by the multiple colored set 
of triangles bounded by green lines. The green lines are the edges to the rarefaction wave, 
also referred to as a rarefaction fan, and the colored triangles represent the different states 
within this fan. Shown on the top, the rarefaction wave has a the left speed of S2 = 0.3972 
and a right speed of = 0.9333. 

On the bottom, there are three separate panels. The left panel displays the Riemann 
invariant space with the three states associated with the Riemann problem under consider- 
ation. The left state uj^, shown as a green dot, is connected to the middle state um, shown 
as a yellow dot, by the 1-shock curve. Also, the middle state is connected to the right state 
Ur, shown as a red dot, by the 2-rarefaction curve. The other two panels show the density 
and velocity profiles after one unit of time. Here one can see, going left to right, the left 
state, the 1-shock wave, the middle state, the 2-rarefaction wave, and the right state in 
these graphs. 

4.2. Time Dilation Between Space Time Cells 

With the capability of solving the Riemann problem for the relativistic compressible 
Euler equations in a unitary frame, like Minkowski spacetime, this section is dedicated to 
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I — Ax — I I — Ax — I I — Ax — I 
Figure 4.5. Effects of time dilation 

solving it in non-unitary frames. A non-unitary frame has the (coordinate) speed of hght 
different from one (i.e. V AB 7^ 1). This factor V~AB, determined by the metric components, 
only changes the speeds of the waves, and when considering solutions to the Riemann 
problem, it has no effect on the states themselves |JGB06] . Since the locally inertial 
Godunov method has a fixed spatial distance for all the frames, this factor determines 
how time is sped up or slowed down relative to the unitary frame. One way to view this 
phenomenon is this factor stretches or shrinks the height of our spacetime cell relative to 
other frames, as illustrated in Figure HT2l This figure shows, from left to right, a frame of 
speed two, a unitary frame, and a frame of speed one-half. One problem we face is one unit 
of time has different meanings in different frames. For example, one unit of time in the left 
frame in Figure 14.21 corresponds to a half a unit of time in the middle frame and a fourth 
of a unit of time for the right frame. In other words, we need a frame to give meaning to 
the quantity of time t, and we choose the unitary frame to be this time keeper for all of our 
simulations. Another problem this factor causes is the need to unify all the times across the 
board, making sure to not break any one frame's Courant-Friedrichs-Levy (CFL) condition 
|LeV92j . To accomplish this unification, the frame where time moves the fastest sets the 
change in time At. This reduction of time in the other frames is represented in Figure 14.21 
as the dashed line. The goal of this section is to determine the effect of shortening the time 
on all the other frames. 

To proceed, we need to distinguish between two different spacetime cells. In the last 
section, we dealt exclusively with a cell containing the whole Riemann problem posed, which 
we will refer to this cell as a Riemann cell. Now, we consider a cell containing an entire 
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Figure 4.6. The Godunov Cell 



Godunov step, and we refer to it as a Godunov cell. In order to implement the Godunov 
step, there has to be at least three states ul, uc, and ur, which poses two Riemann 
problems, and the step is performed on the center state uc', therefore, the Godunov cell 
and the Riemann cell are staggered against one another. Note that if we do not violate the 
CFL condition then the two Riemann problems cannot interact, and the solution u{t,x) is 
completely determined. Define to be the zero speed state for the left Riemann problem 
and the corresponding one for the right Riemann problem. All of these details are 
displayed in Figure 14161 

The following theorem expresses if the time in a Godunov cell is shortened, the resulting 
average is an affine combination of the original average and the center state, based on the 
ratio of the original and new time change. 



Theorem 4.2.1. Let 
(4.42) 



At 



Ax 
2y/AB 



represent the maximum time in a Godunov cell before the CFL condition is violated. Lf the 
change in time is shortened from At to At < At in a Godunov cell containing the solution 
to the Riemann problems u{t,x), the average across that grid cell at time At, u{At), and 
at time At, u{At), are related by 



(4.43) 



u{At) = \u{Ai) + (1 - X)uc 



where X = ^ <1 is the ratio between the two times. 

Proof. Suppose we have a solution to both Riemann problems u{t, x) in the Godunov 
cell, along with a maximum time At (14.42P and a shorter time At < At. Since each 
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half of the Godunov cell contains a disjoint self-similar Riemann problem, we prove the 
claim on one half of the grid cell. Without loss of generality, we study the left half (i.e. 
< X < A.x/2) with the origin located at the bottom left of this cell. We start by 
constructing the relationship between the function u{t, x) at the two times. The fastest 
speed in the Godunov cell (i.e. the speed of light) is referred to as a = V AB. By the 
self-similarity of the Riemann Problem, the following relationship holds 



(4.44) u{At,x) 



for aAt <x<^. 



u{X-^At, X-^x) = u{Ai, X-^x) for < X < aAt 
uc 

Figure HiT] displays these details, where the center state uc is sectioned off by the speed of 
light of the left and right Riemann problems. 

u{Ai, X'^x) 



At 




At 



I Ax 1 

Figure 4.7. The effect of shorten the time step within a Godunov cell 



The average at the time At is directly calculated as 

f-aAt 



Ax 

u{At) =— / u{At,x)dx 



(4.45) 



Ax 
2 

"Ax 







Ax 



u{Ai, X ^x)dx + / ucdx 







Ax 
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At 



Ax 
2 



A u{At,y)dy + uc -aAt 



Ax 



Since ^ = aAt, this implies 



(4.46) 

and ()4.45p becomes 



Ax , \ Ax / aAt\ Ax, 
_-„At =- 1--^ =-(l-A), 



(4.47) 



u{At) = Xu{Ai) + (1 - X)uc, 



proving the claim. 



□ 
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Since the Godunov step only considers the zero speed states and and the factor 
y/AB has no affect on scahng a speed of zero, we expect the affect on the Godunov step 
of this time dilation to be minor. A quick calculation for the left half of the Godunov cell 

shows 

u{At) = Xu{Ai) + (1 - X)uc 



(4.48) 



2Af 

= (1 - \)uc + xuc - ^^[fi^c] - f{ul:)] 

9 At 

= nc-A— [/(txc)-/(n^)] 

A.T 

= uc-^[f{uc)-f{u^)]. 
Hence, the only effect of a reduced time is to shorten the input time in the Godunov step. 



4.3. Locally Inertial Godunov Method 

Equipped with all the necessary tools, we state the algorithm of the locally inertial 
Godunov method. This method is a fractional step scheme started by choosing some pa- 
rameters. In particular, we choose a minimum radius rmin, a maximum radius rmax, the 
number of spatial gridpoints n, and a start time to. In our simulations, the number of 
spatial grid points n is chosen to be a power of two (i.e. n = 2^ for some k). Prom these 
parameters, the mesh width Ax is determined to be 

(4.49) Ax = """"^ " ""'"^'^ 



n - 1 

and is fixed throughout the scheme. Let {xi,tj) represent a mesh point in an unstaggered 
grid defined on the domain 

(4.50) D = {rmin <Xi< Tmax, tj > *o}- 

The spatial points are defined as 

(4.51) Xi = rmin + l)Ax for i = 1, . . . , n. 



Unlike the mesh width, the time step or the mesh height. At, changes from one time step 
to the next because there is no way to determine beforehand the smallest At satisfying the 
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Figure 4.8. The Riemann cell K 



CFL condition for every time step. So for every time tj, a new time step is computed by 
(4.52) Atj = min 



I AijBij J 

where the minimum is taken over all the spatial gridpoints at time tj of the metric Aij = 
{Aij, Bij), where these entries are defined shortly. Starting at to, our temporal mesh points 
are defined as 

j 

(4.53) tj = to + XI ^^'^ for j = 1, . . . , oo. 

k=l 

We assume at our current time tj for j > there exists a solution u{tj,x) and A{tj,x) 
for {tj,x) G D. This solution is either provided as the starting solution at to or from the 
last iteration of the locally inertial Godunov scheme constructed inductively. To implement 
the method, this solution is discretized into piecewise constant states. Discretizing the 
conserved quantities u{tj,x), let tiAx be given by piecewise constant states Uij at time 
t = tj' as follows: 

(4.54) UAx{t:x) = Uij = u{tj,Xi) for Xj_ < x < Xi+,t = fj. 

For notational convenience, we denote and Xi- = X- 1 throughout this paper. 
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We define the grid rectangle Rij so the mesh point {xi-,tj) is in the bottom center of 

it, 

(4.55) Rij = <x<Xi, tj <t < tj+i}, 1 < i < n + 1, j > 0, 

which is diagrammed in Figure 14.81 Each grid rectangle is a Riemann cell, containing 
a solution to a distinct Riemann problem. We are limited to solving Riemann problems 
within Riemann cells having a constant speed of light, as discussed in the last section. 
To this end, the metric source A = {A, B) must be approximated by a constant value, 
denoted Ay, in each Riemann cell Rij throughout the simulation. These constant values 
are established by setting 

(4.56) AAx(t, x) = Aij = A{tj,Xi-) for (t, x) e Rij. 

This approximation makes Aai discontinuous along each line x = Xj, i = 1, . . . , n, at each 
time step t = tj. 

To implement the Godunov step, we need boundary profiles at the left and right bound- 
aries along with the initial profiles at time to because the Godunov step is a three point 
method, and the points xi and x„ need left and right partners, respectively, to pose the 
boundary Riemann problems. These boundary profiles are used to implement the boundary 
Riemann problems and are referred to as ghost cells. The left and right ghost cells, located 
at the points xq and respectively, must be consistent with our numerical solution to 

the Einstein equations around these boundaries. Any inconsistency in these boundary con- 
ditions would result in errors propagating into our solution, corrupting the data; therefore, 
data is needed for the left ghost cell uqj and Aqj along with the right ghost cell Un+ij and 
A„+ij that are solutions to the Einstein equations synchronized with the data close to the 
boundary. Figure 14.91 displays the location of these ghost cells. 

The discontinuities of the metric Aax are staggered relative to the approximate solution 
UAx as illustrated in Figure 14.91 This staggering puts constant metric values within each 
Riemann cell and constant conserved quantities states at the bottom of each Godunov 
cell. Constant conserved quantities uax in each Godunov cell and a constant metric Aax 
in each Riemann cell enables us to pose Riemann problems in locally inertial coordinate 
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Figure 4.9. Staggering of the metric A and the solution u 




frames where we are capable of solving the relativistic compressible Euler equations. More 
specifically, there is a Riemann problem at the bottom center of each Riemann cell Rij 



(^•^'^) I UL = Ui-l,j X < Xi- 

UR = UiJ X > Xi- 

Let uf^{t,x) denote the solution of (I4.57P within the Riemann cell Rij, and define 

(4.58) u^x{t, x) = u^j^{t, x) for {t, x) G Rij 

as the Riemann problem step of the fractional step scheme. 

Equipped with the solutions to the Riemann problems in each Riemann cell Rij, we 
implement the Godnunov step to obtain the average of fiuid variables across the in- 
tervals [xj-_,Xj+] at the next time step tj+i. Since the metric A is different on both sides 
of Xi, separate averages must be taken over the left and right half cells and combined to 



ifj and u^j 



obtain the true average. In particular, let u^^ and ti^ be the average on the left and right 



half cells, respectively. Also, let ul^ = u^j {tj^^, Xi-) and = represent 
the zero speed states left and right Riemann problem, respectively, as shown in Figure 14.61 
To perform the Godunov step on the left half cell, we compute 

(4.59) ufj = Uij - —{f{Aij,Uij) - f{Aij,u^)}, 
and do the same for the right half cell 

r> 2At r, 

(4.60) u^j = u,j - — {/(Ai„ nf ) - f{Mj,Uij)], 
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where the 2 accounts for the half cell calculations. Taking the average of these results leads 
to 

(4.61) nj = ^{ufj+ufj}, 

defining our Godunov step of the method. 

We proceed to define the ODE step. Let u{t,uo) denote the solution to the following 
initial value problem 

ut = G{Aij,u,x) = g{Aij,u,x) - A' ■ V xf{Aij,u,x), 

(4.62) 

n(0) = uq, 

where G{A,u,x) = {G^,G^) takes the form 

,4.63) G» = -IVAB (^) {2(i + 1) - - . 



(4.64) G^ = --VAB j ^ |4.^ + (- - l)(c^ + v^) + -(a^ - v^ypx^j . 

We define the approximate solution UAx{t,x) and AAx{t,x) analytically to derive the 
piecewise formulas used to update the numerical scheme and to be used in the convergence 
proof of the next chapter. The conserved quantities are defined by the formula 

(4.65) UAx{t,x) = ui^{t,x) + f {G{Aij,u{C-tj,u^^{t,x),x)}d^ 

Jtj 

Therefore, UAx{t,x) is equal to u^^{t,x), the solution to the Riemann problems, plus a 
correction term from the ODE step of the method. The metric is derived from the definition 
of the mass 

(4.66) Max{x, t) = Mr^,^ + ^ / u%,{r, tydr. 
In terms of these equations, define the metric as 

(4.67) A,A^,t) = l-'-M^, 
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and 



(4.68) i?Ax(x,t) = Br, exp f |MaxM)^^J ^ rl,^(nA.(r, dr. 

Finally, in order to update the metric and conserved quantities, we use the Riemann 
problem averages Uij to replace the Riemann problem solution u^^{t,x) and perform nu- 
merical integration on the analytical equations (|4.65p - (|4.68p . This process leads us to define 

(4.69) ^Xij+i = Uij + |G(^(Aij + Ai+ij),u{(_ - tj,Uij,x))^ Atj. 
The mass is 

(4.70) M,,,+i = M,,_ + ^ ^ {u%Jxk-,t,+i)xl_Ax) , 



2 

k<i 



with 



(4-71) u'i^{xk-,tj+i) = ^{u^_ij+i + u^j+i}, 
and the metric becomes 

(4.72) A,^, = 1 - 

Xi— 

and 

(4.73) Bij+, = i?,_e^ 
where 

(4.74) r = - 1 + ^Ti|(nA.(x,_,t,+i))Ax| , 



. k<i 



with 



(4.75) UAx{xk-,tj+i) = ^{uk-i,j+i +Uk,j+i}- 

Note that since the metric is staggered relative to the conserved quantities, we use the in 
between values, like and u/^x{xk-^tj-^-l) in the update step. Let Ajj+i = (Ajj+i, i?jj+i) 
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denote the constant value for Aa^- on This concludes the update step and completes 

the definition of the approximate solution uax and Aax by induction. 

To summarize the method, after the setup of the gridpoints, Riemann cells, initial pro- 
files, and ghost cells, the locally inertial Godunov method constructs the solution inductively 
with four major steps: a Riemann problem step, a Godunov step (with time dilation), an 
ODE step, and an update step. The Riemann problem step is described in equations (j4.57p 
and (j4.58p . Formulas (j4.59p - (|4.6ip denote the Godunov step. The ODE step is detailed in 
(ji^ - (j05]) . Finally, equations (|ir69]) - (|T75|) express the update step. 
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CHAPTER 5 

Convergence of the Method 

The focus of this chapter is proving the main theorem of this thesis. This theorem states 
that if a solution (uax, Aax) — >■ {u, A) using the locahy inertial Godunov scheme converges 
and has a total variation bound at each time step, then (u, A) are weak solutions to 

ut + /(A,n)^ = ^(A,^,^), 

(5.1) 

A' = h(h, n, x), 

which is weakly equivalent to the Einstein equations ()2.9p - ()2.12p . This proof is a modifica- 
tion of the Groah and Temple argument using the locally inertial Glimm scheme |JGB06] , 
with a few differences. The main difference is the solution update at each new time step. 
In this paper, an average of the fluid variables is taken verses a random sampling, but the 
steps leading up to this update step are the same. Another difference is the assumed total 
variation bound is used to bound the Riemann problem solutions, as opposed to Groah and 
Temple used wave strengths to bound the Riemann problem solutions. Also, the time steps 
are now variable instead of constant. The last difference is the inclusion of right boundary 
data along with the left boundary data because of the limited extent in space of a computer 
simulation. 

There are two main points to the proof. The first point is to show the discontinuities 
in the metric A along the boundary of Riemann cells are accounted for by the inclusion of 
the term 

(5.2) A' ■VAfiAij,u,x) 

in the ODE step (|4.62p . The second point is to show the jump in the approximate solution 
UAx along the time steps are of order Ax. In their work |JGB06] . Groah and Temple did 
not need the convergence and total variation assumptions because with the Glimm scheme, 
these assumptions are proven as long as there exists a total variation bound on the initial 
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data, a truly remarkable feature of the scheme. In this thesis, these assumptions are shown 
numericahy by the simulation results in the subsequent chapters. Thus, our theorem is 
perfectly suited to our numerical simulation: we numerically establish convergence and a 
total variation bound, and the theorem of this section proves that assuming these conditions, 
we can conclude convergence to a weak solution of the Einstein equations. 

5.1. Convergence to a Weak Solution 

The main theorem of this thesis is the following 

Theorem 5.1.1. Let u/^xit,x) and A^^it, x) be the approximate solution generated by 
the locally inertial Godunov method starting from the initial data UAxito^x) and Aj\x{to,x) 
for to > 0. Assume these approximate solutions exist up to some time tend > to o,nd converge 
to a solution (uax, Aax) — ^ iu,A) as Ax —J- along with a total variation bound at each 
time step tj 

(5-3) T.y.[,_,,_]^A.(t,v)}<^, 

where T.V.[r^.^^^.^^_^j{u/\xitj , ■)} represents the total variation of the function UAxitj^x) on 
the interval [rmim^max]- Assume the total variation is independent of the time step tj and 
the mesh length Ax. Then the solution {u,A) is a weak solution to the Einstein equations 

Proof. Suppose we have approximate solutions {u^^, A^^) obtained by the locally in- 
ertial Godunov method that satisfy the hypothesis of the theorem. Having a total variation 
bound at each time tj places a total variation bound on the inputs to all the Riemann 
problems posed at that time. In |JGB06] . Groah and Temple show a total variation bound 
on the inputs implies a total variation bound on the solution to the Riemann problem for 
any time t such that tj < t < tj+i. By the self similarity of the solution to the Riemann 
problem, this result also implies a total variation bound for any space coordinate within 
the Riemann cell. More specifically, we have the following bounds: 



(5.4) 



T.V.[oc,_^,oc,]{uAx{t,-)} < V, 
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and 

(5.5) r.K[i^,,^.^,)^A.(-,x)} <y, 

for any x and t within the Riemann cell Rij. 

All the functions /, G, and g derived in [JGB06] are smooth, and it is the metric that 
is only Lipschitz continuous. The smoothness of these functions is used throughout this 
proof. 

Let T = tend — to be the overall time of the solution, and for each mesh length Ax define 
the minimum time length 

(5.6) At = mm{Atj} 

3 

as the minimum over all the time lengths defined by (|4.52p . By definition, this time length 
is proportional to the mesh length. At oc Ax, implying 0{At) = 0{Ax), and there exists a 
constant C bounding all the time lengths, Atj < CAt for all j. Throughout this chapter, let 
C be a generic constant only depending on the bounds for the solution [to, t(.nd\ x [rmin-, fmax]- 
This variable is created to unify all the time steps, and more importantly, used to calculate 
the maximum number of time steps needed to go from to to tend- 

We now follow the development of Groah and Temple in |JGB06] . Recall, u^{t,x) 
denotes the collection of the exact solutions in all the Riemann cells Rij for the Riemann 
problem of the homogenous system 

(5.7) ut + /(Aij,u)^ = 0. 

So u^{t,x) satisfies the weak form of this conservation law in each Riemann cell 

0= {-u^^ipt- f{Aij,u^^)ip^}dxdt 

■J Jn^j 

+ j {u^^{tj+ux)ip{tj+i,x) -u^^{t1. ,x)Lp{tj,x)\dx 

(5.8) ■'^^ 

+ / (t,Xi))(/?(t,x,) 
-f{Aij,u^^{t, Xi-i))ip{t, Xi_i)} dt, 
where 93 is a smooth test function with Supp{ip) C [to, ^end) x [o-i b] for a < Vmin < T~max < 
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Remember, u(t, uq) denotes the solution to the ODE 

ut = G{Aij,u,x) = g{Aij,u,x) - A' • V Af{Mj,u,x), 

(5.9) 

■u(O) = Uq. 

Therefore, 

(5.10) u(t,UQ) = UQ+ I {fif(Aij-,n(^,no),x) - A' • VA/(Aij,n(^,no),x)} d^. 

JO 

Also, recall uax denotes the approximate solution obtained using the fractional step method. 
Since our fractional method takes the Riemann problem solution and feeds it into the ODE 
step, u^x is defined on every Riemann cell Rij as 



(5.11) 



ft 

u/s,x{t,x) = u^'^{t,x) + I {g{Aij,u{(, - tj,u^^{t,x)),x) 

Q f 



This expression implies the error between the approximate solution and the Riemann prob- 
lem solution is on the order of Ax; a fact that is repeatedly used throughout the proof. 

Define the residual e = £{uax, Aax, f) of uax and Aax as the error of the solution in 
satisfying the weak form of the conservation law (|5.ip by 



(5.12) 



:iuAx,AAx,v)= / {-UAx'Pt - fiAAx,UAx)Vx - giAAx,UAx,x)(p} dxdt 



h-h 

=n+l 



/ {-UAx<^t - f{Aij,UAx)Vx - g{Aij,UAx,x)ip} dxdt 



-h- h, 



where 



(5.13) UAx{to,x)dx = '^ UAx{t'^,x)dx, 
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and 



(5.14) 



to 



Yj / {•/'(^«i'^Aa'(t,r-+,„))¥'(*,'^™„) - /(Aij,nAx(i,r+„x))(^(t,r+„^)} dt, 



The expression Yli=i~j^ denotes a double sum where the index i runs across ah the spatial 
gridpoints, and the index j runs across all the temporal gridpoints. Remember, n is the 
number of spatial gridpoints, and there are n + 1 Riemann cells as depicted in Figure 
14.91 Our goal is to show e{uAx,-AAx,^) = 0(Ax) because if the approximation converges 
(uAx) Aax) — ^ {u,A) as Ax — )• 0, then the limit function satisfies the condition of being a 
weak solution to the Einstein equations e{u,A,Lp) = 0. 
Substituting (j5.1ip into (|5.12p gives us 

i=n+l „ „ 

{-U^xft - fiAij,UAx)Vx - giAij,UAx,x)ip 

^^■^^^ -(ft [ [g{Aij,u{^ - tj,u^^{t,x)),x) 

- (Ajj, n(^ - tj,UAxit, x))) ■ A'ax] d^} dxdt - h - h- 



Define 



(5.16) 



Iij{t,x)= I [9{Mj,H^ -tj,u^^{t,x)),x) 
d f 

- {Aij,uiC - tj,u^^it,x))) ■ A'^ d^ 
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Plugging the weak form of the conservation law (j5.8p of each grid rectangle into (j5.15p gives 
us 



(5.17) 



i=n-\-l 



X] / / Wx[fiAij,U^^) - f{Aij,UAx)]- giAij,UAx,x)ip 

—(ftl}j{t, x)} dxdt 
i=i,j 

i=n+l „ 

h- ^ {fiAij,U^^{t,X^))ip{t,Xi) - f{Aij,U^^{t,X^-l))Lp{t,X^-l)} dt. 



i=n+l 



Note 



(5.18) 



which implies 



(5.19) 



i=n+l 



^ f[f{Mj,u^^) - f{Aij,UAx)]dxdt 

— 1 o Rii 



< C At' Ax (1^) (n + 1) = 0{Ax) 

where the number of time steps is proportional to T/ At and the number of space steps is 
0{l/Ax) by (ICTIl . 

Since u^{t'^^x) = UAx{t'j ,x), the following sum is rearranged to become 



i=n+l 



(5.20) 



-h- / {'^Axitj+i,x)ipitj+i,x) -u^^{t^,x)ip{tj,x)^ 
i=i,j ■'^^ 

|wAx(ij',2;) - UAx{t],x)^ ip{tj,x)dx 
(p{tj,x) luAx{tj',x) - UAxitJ,x)j dx 

j-ty - min 

+ X] / 'Pitjix)\^UAx{tJ,x) -UAx{tJ,x)jdx, 



dx 
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where the term UA.x{tj,x) is added and subtracted to isolate the jump in the solution uai 
across the time step tj. We define this jump ei = £i{ua.x, Aax, v) a-s 

ip{tj,x) ^UAx{t^,x) - UAx{tJ,x)j dx, 
and this definition allows us to rewrite (|5.17p as 

i=n+l „ „ 

e=0(Ax) + ei+ ^ {-g{Aij,UAx,x)(p - (ptllj{t,x)} dxdt 



ip{t,x) }^UAx{t-,x) -u^^{t-,x)j 

min 



dx 



i=n+l 



i=l,j •'Rj 

But the last sum is rearranged to cancel the boundary conditions as follows: 

i=n+l 



~^2- V / {fiAij,u^^it,Xi))^{t,Xi)-f{Aij,u^^{t,Xi-i))ipit,Xi-i)}dt 

i=n „ 

= X] / {fi^i+iJ,UAxit,Xi)) - fiAij,u^^{t,Xi))} ^p{t,Xi)dt 
(5.23) ^=1 J ' 

+ / {f{Mj,u^xit,xo))-f{Aij,UAxit,xo))}^p{t,xo)dt 

~^yZ {/(A„+l.j,uf^(t,a;„+i)) -f{An+l,j,UAx{t,Xn+l))}^pit,Xn+l)dt, 



where 



(5.24) 



/ {f{M,j,u^{t,XQ)) - f{Aij,UAx{t,XQ))] (p{t,xo)dt 



<||^ILCAt2(|-^ =0(Ax) 



and similarly 
(5.25) 

/ {f{-^n+l,i-,U^{t,Xn+l))-f{An+lJ,UAx{t,Xn+l))}^{t,Xn+l)dt 



0(Ax) 



Note that the resulting double sum in ()5.23p lost a term, resulting in only n terms. 
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To simplify the I^j term, we add and subtract a term deviating from it by an order of 
Ax, use integration by parts on the new term, and with the result add and subtract another 
term to reduce the expression further. To this end, let 

(5.26) 

i=n+l „ „ „t 

Ias= X] / / {9{^ij,u{C-tj,u^^{C,x)),x)-g{Aij,u{^-tj,u^^{t,x)),x) 
1=1 J •' •'^j •'tr 



iAij,u{^ - t, u^^it, x))) ■ A'^J d^dxdt. 

From the total variation bound on the Riemann problems and the smoothness of /, this 
term is bounded by 

i=n+l .. .. ..t 

\Ias\< E / / ll^*lloo/ C T.V.^,^_^^,^^{ui,:,{;tj)}didxdt 

i=l,j "^^ij 

(5-27) <||^,||^CAt2Ax^r.y.[,^,„,,_]{uA.(-,t,-)} 

j 

<CF||^i||^AxAi2|- = 0(Aa;2), 
and the above procedure reduces the term to 
(5.28) 

~ ^tlljit,x)dxdt = Ias- H / / {5(Aij,M(C-ij,wf^(C,a;)),a;) 

J ^ Rij i=l j ^ Rij ^tj 

- -qJ^ (Aij, n(e - ij, ?iAx (C> x))) ■ A'^^} dxdt 

i=n+l f tj+i 

= 0(Ax2)- ^ / lip{tj+i,x) [g{A,^,ui^-tj,u^^{^,x)),x) 

- 1^ (A,„ u(C - t,,uZ (e, x))) ■ A!^,] di 

- j ^[g{^ij,UAx, x) - -^iAij,UAx) ■ ^Ax¥^ I dx 

= 0{Ax'^)- }2 / \v>itj+i,x) [giA^^,ui^-tj,ul^itj+i,x)),x) 



OA 



{Aij,u{^ - tj,u^^{tj+i,x))) ■ A'^^] dt + h + h, 
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where 
(5.29) 



i=n+l 



g{Aij,u{i-tj,u^^{i,x)),x) - —{Aij,u{i-tj,u^'^{tj+i,x))) • A'^^ 



+ ^ {Mj,u{i - tj,ul^{C, x))) ■ A'^J dej dx, 



and 
(5.30) 



i=n+l 



Ri 



g{Aij,UAx,x) - -g^iAij,UAx) ■ A'^^ 



dxdt. 



Again by smoothness and the total variation bound, we have 

i=n+l 



\h\ < llv'lloo X] ^ ^-^-[x,-!,!:,] {^^Aa;(-,ij)} AxAt 



(5.31) 



< ||<^||^CAxAt J^r.F.[,_,,_] {uU;h)} = M\^CV\xAt— = 0(Ax). 



Substituting (I5.23P and (I5.28P into (I5.22P along with using ()5.1ip as an identity leaves us 
with 



(5.32) 



i=n+l 

£ = 0(Ax) + ei - ^ 



Qf 

V9— (Aij,UAx) • A'^^dxdt 



i=n „ 

/ ip{t,Xi){f{Ai+ij,u^^{t,Xi)) - f{Aij,u^^{t,Xi))}dt 



1=1,3'' ^1 



The second sum represents the jump in the flux function /, resulting from the discontinuities 
in the metric A, and the first sum is the addition to the ODE step (|5.9p specifically designed 
to cancel these jumps in the flux. 

To see how the cancelation works, we perform a Taylor expansion on the test function, 
and we add and subtract terms deviating by order Ax. The first sum in (|5.32p is expanded 
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as 



(5.33) 

j=n+l 



^ / / if{xi,t) — {Aij,UAx) ■ -^Axdxdt + 0{Ax) 

= ^ I '^(^i' / \ ^(Aij, MAx) • A'^^. - — (Ajj, u^^) • A'a^ \ dxdt 

+ E "Pi^^^i) {|^(A^,,n^^)-A'^,-|^(Ai„u^i:(xi,t))-A'^,|dxdi 

1=71+1 „ ,-Xi 

+ / ip{Xi,t) —{AAx{x + —,tj),U^^{Xi,t))-A'^Jxdt + OiAx) 



From the smoothness of /, each of the first three sums in equation (I5.33P are 0(Ax) for 
the following reasons: tlie first sum is order Ax from the ODE step in the definition of the 
approximate solution u^x (I5.11P . the second sum is order Ax'^ by the total variation bound 
on solutions to the Riemann problems, and the third sum is order Ax by the Lipschitz 
continuity of the metric A. After these bounds are established, (j5.33p reduces to 

(5.34) 

df 



^ V'^(Aij,nAx) • A'^^dxdt 

I ^{xi,t) f ^{AAx{x + ^,tj),u^^{xi,t))-A'^Jxdt + 0{Ax) 

— 1 j ^i— 1 

=!^^ f f^^ df Ax 

V / f{xi,t) j^{Aax{x + —,tj),u^^{xi,t))dxdt + 0{Ax) 
tl^j JRj Jxi--, ox I 

/ ^{t,X,){f{A,+lJ,ul^{t,X^+))- f{Aij,U^^{t,X^+))}dt + 0{Ax). 



i=l,j •'^3 
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Plugging this result ()5.34p into (j5.32p gives us 

e = 0(Ax) + £i 

(5.35) r 

- X] / <^(*>a;„+l) {f{An+2,j,U^{t,Xn+l)) - f[An+l,j,U^{t,Xn+l))] dt, 

where one term remains due to the mismatch in the number of terms in the spatial sum. 
Clearly, this last term is 0(Ax). 
So the residual boils down to 

(5.36) e('UAa;, Aacc, 'f) = ei{uAx, Aax, ^) + 0(Ax), 
with all that remains to show is 

ip{tj,x) |'UAx(i/,a;) - u/^x{t~,x)\ dx = 0{Ax). 

■_ . . ■ m. 7 n. 



To estimate Si , we break up the sum by each time step tj and define 



(5.38) 



ip{tj,x) \uAx{t~j,x) - UAx{t~,x)j dx 

n 

j (pitj,x)^UAx{t'j',x) -UAxitJ,x)^ dx, 



with Xi+ = X- 1 1 and Xi_ = x- i. 

Recall, the approximate solution for the new time step is computed by the Godunov 
step, using averages at the top of each Riemann cell Rij. In particular, the solution at each 
new time step is 

(5.39) UAx{tj',x) = u{tj - n(tj), x)) 
where 

(5.40) u{tj) = ^ ul^{tj,x)dx 

To finish the proof, a lemma is needed, which is proven at the end of this chapter. This 
lemma states the difference of the ODE step taken on an average verses the solution to the 
Riemann problem across the top of the Riemann cell is bounded by the total variation of 
the Riemann problem. 
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Lemma 5.1.1. Let represent the solution of the Riemann problem in the Riemann 
cell Rij-i and UAx{t) denote the average of the Riemann problem solution across Riemann 
cell. Let u be the solution obtained by the ODE step ^5.9\) and if be a smooth test function. 
Then the following bound holds 



(5.41) 



{u{tj — tj-i,UAxitj), x) — u{tj — tj-i,u^^{tj,x),x)^ (p{tj,x)dx 



for some constant C. 

Using Lemma [5.1.11 (j5.38p is rewritten as solutions to the ODE step (15. 9p and bounded 

by 

(5.42) 

£{ = 22 fitj^x) - tj^i,u{tj),x) - u{tj - tj^i,u^^{tj,x),x)} dx 

< ciiv^lL AxAt^r.y.[,^_,,^^]{n^^(-,t,)} = c\M^AxM r.y.[,_,,_]Ki:(.,t,)} 

i 

By the total variation bound on U/\x{tj, the residual is bounded by 

(5.43) ei < ^CIIv^lL AxAt T. y. n^^ (•, t,-)} < C^^AxAtV = 0{Ax). 

Therefore, e = 0{Ax) and the proof is complete. □ 

To prove Lemma [5.1.H a preliminary result is needed: given a function on a set of points 
the difference of the function between any point and the average is bounded by the total 
variation of that function on the set. This result is provided by the following 

Lemma 5.1.2. Let u(x) be a function on the set and 

(5.44) u = — j u{x)dx 
be the average of u on this set. Then we have 

(5.45) \u - u{x)\ < sup \u{xi) - u{x2)\ < T.V.i^^_^^^^]{u{-)}. 

Xl,X2'^[Xi,Xi+l] 
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Proof. The second inequality is true by the definition of the total variation 



(5.46) 



sup \u{xi) - U{X2)\ < T.V.[^^_^^^^]{u{-)}. 

xi,x2e[xi,Xi+i] 



To prove the first inequality, we assume it is false to obtain a contradiction, so suppose 
there exists x* G such that 



(5.47) 



sup \u{xi) — U(X2)\ < \u — u{x^:)\. 

xi,x2e[xi,Xi+i] 



Relabel the u-coordinates by an isometry (p : u ^ v that maps the point u(x*) to the origin 
in the v-coordinates (i.e. ^p{u{x^)) = 0), and the vector u — u{x^) in the direction of the 1st 
coordinate v^, as show in Figure [5Tl 



U2 



V2 




d 



^(u) = V 



O = $(n(x*)) 
Figure 5.1. The isometry ^ : u ^ v 



Vl 



Since the average of a collection of points is independent of the coordinate system in 
which they are labeled in, we have 

\ r^i+ \ r^i+ ( \ r^i+ \ 

(5.48) V = —— / v{x)dx = —— / (p{u{x))dx = ip —— / u{x)dx = <p{u) 

^^Jxi^ ^^Jxi- Jxi- I 

The following inequality holds by transforming equation ()5.47p over to v-coordinates 



(5.49) 



\v{x)\ = \u{x) — ti(x*)| <\u — n(x*)| = \v\ Vx G [xi_, Xj+] 



which implies 
(5.50) \i 



Ax 



v{x)dx 



Xi 



< - — / \v(x)\dx < - — / \v\dx = \v\. 
~ Axy,^ ' ' Ax L ' ' ' ' 



This inequality \v\ < \v\ is an obvious contradiction, proving the first inequality in (|5.45p . 

□ 
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Now we prove the lemma used in the proof of Theorem 15.1.11 
Proof of Lemma 15.1.11 Recall the solution to the ODE step has the form: 



(5.51) u{tj -tj-i,u^^{tj,x),x) = u^^{t,x) + / G{Aij,u^^{t,x),x)dt. 

This solution implies the LHS of (|5.4ip is written out as 

{u{tj - tj-i,UAxitj),x) - u{tj - tj.i,u^^{tj,x),x)] ip{tj,x)dx 

{{uAxitj) -u^^{tj,x)) 
{G{Aij,UAxit),x) - G{Aij,u^^{t,x),x))dt] ip{tj,x)dx 



+ 



(5.52) 



+ 



{uAxitj) - UAxitj^x)} >fitj,Xi)dx 



2\ 



{G{Aij,UAx{tj),x) - G{Aij,UAx{tj,x),x)} dt ip{tj,Xi)dx 



+ 0{Ax' 

where the test function in the first term is approximated by a Taylor expansion. By the 
definition of the average function u, the first term is zero. By the smoothness of G, the 
bound (j5.4ip is proven by 



(5.53) 



{u{tj - tj-l,UAx{tj),x) - U{tj - tj-l,U^^{tj,x),x)} ip{tj,x)dx 



< G\\ip\\^AxAt sup {\uAxitj) - UAxitj, x)\} 



< llV'lloo T.V.[^^_^^^^]{uAx{tj,-)}, 

where Lemma 15.1.21 is used to bound the difference between the average and the solution 
to the Riemann problem. 

□ 
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CHAPTER 6 

Continuous Models 

Before simulating the shock wave models, we briefly explore simulating more simple ones, 
which we denote as the continuous models. The continuous models are the pure FRW-1, 
FRW-2, and TOV metrics developed in Chapter [3l and these models are not as complex 
as the discontinuous shock wave models. Our purpose in simulating these models is to use 
exact formulas to test the numerical convergence and the accuracy of our locally inertial 
Godunov scheme. Each of the three models embodies a different time dilation scenario, 
each of which occur in the simulation of the shock wave models. The FRW-1 metric is a 
model in which the (coordinate) speed of light is uniformly equal to one independent of 
time and space; therefore, time dilation and synchronization does not occur in this model. 
In the TOV model the (coordinate) speed of light increases from one side of the simulated 
region to the other, so there exists time dilation between frames in this model. Since the 
TOV metric is static, the synchronization of the clocks has no effect on the construction 
of the solution. In the FRW-2 model, we choose an appropriate integrating factor to force 
the (coordinate) speed of light to be equal to one initially, but in this model the speed of 
light increases uniformly across the entire model as time progresses. This model provides 
the scenario of a dynamical (coordinate) speed of light, so synchronizing time against the 
unitary frame must be handled correctly to obtain numerical convergence. As a side note, 
it was serendipitous that we came across the FRW-2 model; this model allowed us to 
perfect our clock synchronization by providing an ideal test case. Besides running these 
test cases, the simulation of the continuous models are important for correctly handling 
the non-interaction region and ghost cells for the shock wave model. As discussed in more 
detail in the next chapter, emanating from the initial discontinuity is an interaction region 
surrounded by non-interaction regions on each side. These non-interaction regions are the 
continuous models and their numerical accuracy and consistency is paramount to the success 
in simulating the shock wave model. 
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The models we consider in this paper have infinite extent, with the radii going from zero 
to infinity, and we refer to this infinite region of space as the universe of our simulation. 
Even though theoretically space can extend out to infinity, the computer can only simulate 
a finite region of spacetime in the model, so we need to demarcate the minimum radius fmin 
and the maximum radius fmax of the simulated region, along with the number of gridpoints 
n to simulate it. Also, we need to decide when to begin the simulation of our models, so a 
start time to is chosen. We choose to set these parameters as 

(6.1) fmin = 3, fmax = 7, to = 15, n = 2^'^ = 16, 384, 

in all the simulations within this chapter. The units for these parameters are given meaning 
and are discussed in Chapter [SJ There is an extra parameter for each of the models; it is 
the integrating factor constant for the FRW-1 and FRW-2 models and the constant Bq 
for the TOV model. This parameter is a time scale factor that affects the rate at which the 
clocks run in the model, and it is assigned a value as we consider each model separately. 

This chapter is organized into four sections. The first three sections. Sections I6.m6.31 
are dedicated to simulating the FRW-1, TOV, and FRW-2 metrics, respectively. These 
sections start out by reiterating the equations developed in Chapter [3] that are needed to 
run the simulation. More precisely, these equations are used to build the initial profiles and 
ghost cells. With this data set, we use our locally inertial Godunov scheme to simulate these 
models and provide a glimpse of them with various pictures. In Section [6.41 we discuss how 
the errors and convergence rates are computed, and we record the numerical convergence 
in the simulation of all three continuous models. 

6.1. Simulating the FRW-1 Model 

Recall, using the coordinate transformation 



r 



--Vtr, 



(6.2) _ . -2. 

^=|^+4^^r=^+4' 

the FRW metric in standard Schwarzschild coordinates (the FRW-1 form) is 

(6.3) ds'^ = ^-^dP + —^df"^ + f'^d^f, 

1 — t;^ 1 — 
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Figure 6.1. Initial profiles for FRW-1 model 



with the metric components 

1 



(6.4) A{0 = l-v{0\ B{0 



and the fluid variables 

(0.5) = = 
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where the variable ^ is defined as ^ = r/t and is a function of v by the following relation 

2v 



(6.6) 



The integrating factor constant is suppressed, which is equivalent to setting ^'q = 1. 
With this parameter set, the initial profiles and ghost cells for the simulation are constructed 
by equations (j6.3p - (|6.6p along with the standard parameters (j6.ip . The resulting profiles are 
illustrated in Figure [6Tl The graphs for the density, velocity, metric component A, metric 
component B, and the mass function are sectioned into five panels displayed from top to 
bottom and color coded by the colors red, blue, green, yellow, and brown, respectively. 
Each graph is scaled and shifted accordingly in the panel to ensure the top represents the 
maximum value and the bottom represents the minimum value of the graph in the region 
of space under consideration. The numeric values for the minimum and maximum points 
on the graph are shown on the left side. For example, the velocity profile is the blue graph 
second from the top with a minimum value of 0.1010 and maximum value of 0.2476 between 
the region of space [rmimfmax]- For this particular example, the minimum value occurs at 
the left most point r^m) while the maximum occurs at r^ax- The vertical magenta line 
is a marker the user controls to examine a particular gridpoint in the simulation, with the 
corresponding numerical values displayed in magenta on the left side. Another piece of 
information associated with this marker is the scale factor or (coordinate) speed of light 
(V AB) which is shown under the graphs and labeled "SCALE". Looking at the graphs 
in Figure 16.11 all of the profiles are increasing functions except the metric component A 
is decreasing. The graph of metric component A in Figure 16.11 is the mirror image of the 
graph B, which confirms they are recipricals of one another (|6.4p . 

The simulation is run for one unit of time, tend = + 1) with the results depicted in 
Figure [6^21 This figure features the evolution of the fiuid variables (p, v), showing snapshots 
of the density and velocity profiles at the times, going from the left to the right frame, 
t = to, t = to + 0.2, t = to + 0.4, t = to + 0.6, t = to + 0.8, and t = tend- The fiuid 
variables decrease as time progresses, with the density decreasing faster than the velocity. 
Since a positive velocity indicates matter is moving outward, we expect this decrease in 
density corresponding to a decrease in the mass function. The metric has little change, and 
moreover, the solution as a whole has a very similar shape from the initial profiles seen 
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Figure 6.2. Evolution of the FRW-1 model during a unit of time 

in Figure 16.11 Notice the (coordinate) speed of light is identically equal to one, which is 
confirming the fact the metric components A and B are recipricals ()6.4p at all times. 



6.2. Simulating the TOV Model 



Remember, the TOV metric takes the form 



(6.7) 



1 



2g M(f) 
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shock Wave Cosmology Simulator 
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Figure 6.3. Initial profiles for TOV model 



with the metric components 



(6.^ 



A{r) = 1 - 87rg7, B{f) = Bo{r)^ , 



and the fluid variables 



(6.9) 



p{r) 



7 



v{r) = 0, 
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Figure 6.4. Evolution of the TOV model during a unit of time 

where the parameter 7 is a constant dependent on a 

1 



(6.10) 



7 



a 



27rg V 1 + 6cr + 0-2 / ' 



Setting the extra parameter Bq = 1 and using the standard parameters (j6.ip , the initial 
profiles and ghost cells are built from these equations (j6.7p - (j6.10p (with a = 1/3), and the 
resulting profiles are shown in Figure 16. 3[ Notice the ghost cells can be set initially and do 
not need to be changed since the TOV metric is static, independent of time. This model 
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contrasts with the evolutionary nature of the FRW-1 model in Figure 16.11 In the TOV 
model, the density is a decreasing function, while the velocity and metric component A 
are both constant functions. The metric component B and the mass function are always 
increasing functions since they are computed as the integral of positive values as seen by 
equations (j4.68p and (j4.66p . Since the metric A is constant and the metric B is increasing, 
the (coordinate) speed of light is increasing across the universe, so time dilation between 
the different frames occur in this niodGl. In particulcir, the spGGd of light is 0.58 set r-m.^n 

and 

0.89 at rmax- 

Running the simulation for one unit of time, the evolution of the TOV metric along 
this time frame is shown in Figure 16.41 This metric is static, so we expect it to remain 
unchanged throughout this simulation. The density, metric B, and the mass profiles are 
unchanged. The constant profiles of the velocity and metric A appear to have changed. This 
appearance is due to the small numerical errors and the auto-zoom feature of the graphing 
tool which magnifies these errors. 



6.3. Simulating the FRW-2 Model 



Remember, using the coordinate transformation 



r =Vi', 



;r, 





the FRW metric in standard Schwarzschild coordinates (the FRW-2 form) is 



(6.12) 



1 



1 





1 — V 



with the metric components 






the integrating factor 
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Figure 6.5. Initial profiles for FRW-2 model 

and the fluid variables 

(6.15) p(tf) = ^^, v(t,r) = - = —. 

Unlike FRW-1, the FRW-2 metric relies on the FRW time coordinate t, which is the following 
function of the new variables 
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Figure 6.6. Evolution of the FRW-2 metric during a unit of time 

To match the FRW-1 metric more closely, we choose the integrating factor constant ^'o 
so that the (coordinate) speed of light V AB is equal to one. To this end, we set 



.17) 



This relationship is explained in more detail in the next chapter and is a rearrangement 
of equation ()7.49p . The initial profiles and ghost cells for the simulation are created by 
equations (I6.12p - (16.16p . using the standard parameters (16. ip and ^'o set in (16.17p . These 
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profiles are pictured in Figure 16. 5[ In this figure there is more curvature in the graphs 
of FRW-2 than the FRW-1 model (Figure iG.ip . This curvature produces higher values in 
all the increasing graphs (density, velocity, metric B, and mass) and lower values in the 
decreasing graphs (metric A). 

Allowing the simulation to run for one unit of time, the evolution of the FRW-2 metric 
is shown in Figure 16. 6[ Again, the FRW-2 metric follows the same overall pattern as the 
FRW-1 metric (Figure [6.2p . and the fluid variables are decreasing as time progresses. We 
expect these similarities since both metrics, FRW-1 and FRW-2, are derived from the same 
FRW metric. In the FRW-2 universe, the (coordinate) speed of light changes uniformly 
across the simulated region from 1 at to to 1.0667 at tend- Thus, the FRW-2 model has 
a dynamic (coordinate) speed of light. In contrast to the FRW-1 and TOV models, our 
simulation of the FRW-2 model tests the locally inertial Godunov method in a setting in 
which the (coordinate) speed of light is dynamic. 

6.4. Convergence Results 

This section discusses the numerical convergence results for all three of the continuous 
models. These results are recorded for the FRW-1, TOV, and FRW-2 metrics in Table 
16.11 16. 2( and 16.31 respectively. The tables are organized to show, from left to right, the 
number of gridpoints used, the density results, the velocity results, the metric A results, 
and the metric B results. The variable results are partitioned into two values: the one- 
norm error and the convergence rate. Since we possess exact formulas for the solution of 
all three models, the 1-norm error is numerically computed between the numerical solution 
and this exact solution. The 1-norm is a natural norm to use for conservation laws because 
it requires integrating the function, and the weak form of the conservation law gives us 
information about these integrals. The 1-norm is our chosen method for computing the 
error and showing numerical convergence in this paper. The convergence rate is computed 
by taking the log2 of the ratio in successive errors, enabling us to measure the decrease in 
error relative to the increase in the number of gridpoints. For example, a rate of 1 means 
that using twice the number of gridpoints reduces the error by half. A rate less than 1 
means the error is reduced by less than a half, while a rate greater than 1 means the error 
is reduced by more than a half. Looking at the tables, all the errors are decreasing as the 
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number of gridpoints increase with convergence rates around 1. The convergence rate of 
1 is expected because we are implementing a first order method on continuous solutions. 
These results indicate the simulation is producing an accurate numerical representation 
of all three models, giving us the green light to simulate the discontinuous models in the 
following chapters. 
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9.310e-005 
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Table 6.1. Convergence results for the FRW-1 model 
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Rate 
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N/A 


1.526e-003 


N/A 


3.091e-003 


N/A 
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1.340e-007 


0.99 
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0.98 


7.531e-004 


1 
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1 
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1.894e-004 


1 


2048 


8.450e-009 


1 


2.880e-005 


1 
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1 


9.470e-005 
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Table 6.2. Convergence results for the TOV mode 
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Rate 
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N/A 
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1.195e-002 


N/A 
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1.263e-003 
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0.98 
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1.980C-006 
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1 


3.227e-003 


0.91 
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9.890e-007 


1 


3.136e-004 


1 


1.289e-003 


0.98 


1.899e-003 


0.76 
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4.940e-007 


1 


1.568e-004 


1 


6.240e-004 


1 


7.170e-004 


1.4 
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1 
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1 


3.131e-004 


0.99 


3.704e-004 


0.95 
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1 
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1 


1.579e-004 


0.99 
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0.89 
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1 
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1 


8.030e-005 


0.98 


1.181e-004 


0.76 
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1 


9.790e-006 


1 


3.890e-005 


1 


4.470e-005 


1.4 
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CHAPTER 7 

Shock Wave Model 

In this chapter we perform the numerical simulation of a spherically symmetric general 
relativistic outgoing shock wave, together with the secondary incoming reflected wave. This 
simulation demonstrates the secondary reflected wave is also a shock wave. To reiterate, we 
interpret this result as the numerical resolution of the secondary reflected wave associated 
with the Smoller and Temple shock wave model |ST95|. IST03|, [5T04] . That is, in |ST95j . 
Smoller and Temple (Sm/Te) constructed an exact shock wave solution of the Einstein 
equations consisting of an inner FRW spacetime blasting outward into a TOV spacetime, 
such that the interface between them was a true fluid dynamical shock interface. The metrics 
satisfied equations of state p = ap (FRW) and p = ap (TOV), where a and a were constant, 
consistent with an isothermal scenario. Since the outer TOV solution was inverse square 
in the density, the solution could be interpreted as a blast wave propagating outward into 
a static singular isothermal sphere, c.f. |ST95] . To get exact formulas, the Sm/Te model 
assumed different sound speeds (temperatures) ahead and behind the shock wave. This 
determined cr as a function of a, which we interpret here as having the simplifying effect of 
eliminating the secondary reflected wave in the solution, thereby making the construction 
of exact formulas possible. 

To simulate the secondary reflected wave and develop a picture of these solutions, we 

2 

assume both TOV and FRW spacetimes satisfy the same equation of state, which isp = yp, 
for the pure radiation stage of the early universe. To run the simulation, we use an exact 
formula for the FRW spacetime in standard Schwarzschild coordinates first constructed in 
|STj . detailed in Chapter [3l This transformation puts the FRW metric into the "same 
coordinate system" as the TOV metric, enabling us to match these metrics together, with 
the FRW metric on the inside, to form the FRW/TOV matched metric. Since we have two 
transformations of the FRW metric into standard Schwarzschild coordinates, denoted as 
FRW-1 and FRW-2, we have two matched models, FRW-l/TOV and FRW-2/T0V, at our 



disposal to simulate. This matching gives us exact expressions for the initial data (a point 
of continuity between FRW and TOV metric components, such that the density, pressure, 
and velocity suffer a jump discontinutiy) and boundary data. This data is the information 
necessary to run the locally inertial Godunov scheme, developed in Chapter HI to perform 
the numerical simulation of the interaction region between the two spacetime metrics. In 
particular, since light speed is a speed limit for propagation of signals, outside this lightcone 
the solution should remain FRW (on the inside) and TOV (on the outside). We call this 
region contained within the light cone that emanates from the initial discontinuity the 
region of interaction. Besides simulating shock wave formation, we also compute the cone 
of light and cone of sound, the light and sound information that propagates outward from 
the initial discontinuity. Both the cone of light and the cone of sound emanate and expand 
away in both directions from the initial discontinuity. Because light travels faster than 
sound, the cone of sound must be contained in the cone of light. Since we expect from 
|JGB06j that there is no lightlike propagation in spherically symmetric spacetimes, the 
cone of sound should be the true limit of propagation of information; therefore, outside the 
cone of sound we expect the solution should remain FRW (on the inside) and TOV (on the 
outside). Figure [7T] illustrates all these expectations. This prediction is an important result 
we demonstrate in our numerical simulation. In fact, since the numerical simulation involves 
integration of the metric components at each time step, we know of no mathematical proof 
that the weak solutions constructed by the Groah and Temple theorem |JGB06] actually 
have the sound speed as a propagation speed for the information in the solution. 
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Just like the continuous models, we consider the entire infinite region of the FRW/TOV 
model as the universe of our simulation, and we need to define initial parameters to run 
the simulation. Recall, our chosen parameter to the one parameter family of shock wave 
solutions is the radius of the initial discontinuity fo separating the two metrics. We choose 
to set these parameters as 

(7.1) fmin = 3, rmax = 7, Tq = 5, 71 = 2^^ = 16,384, 

in all the simulations in this chapter, unless otherwise stated. The units to these are 
explained in Chapter [9j Also, throughout this chapter, when we discuss mesh points, we 
use the space coordinate x to match the notation of the locally inertial Godunov scheme as 
denoted in Chapter HJ and when we discuss points within our universe, we use the radial 
coordinate r to explicitly express these points as a radial coordinate of the metric. When we 
discuss the derivative of a function in our solution, we compute these derivatives numerically 
using a standard three point method |BF97j . 

In Section 17.11 we introduce material needed for testing the accuracy and certainty of 
computing the correct solution in our numerical simulation. We start by computing the 
(coordinate) speed of light and sound. We move to tracking the borders for both the FRW 
and TOV side of the simulation, and we conclude this section by finding the boundary 
condition for the TOV side. We setup the FRW-l/TOV matched model simulation in 
Section 17. 2^ determining the initial profiles and boundary conditions we can discretize and 
feed into the locally inertial Godunov scheme. We use this setup in Section 17.31 to run the 
simulation and obtain results. We record the convergence of the entire simulation, using 
the successive mesh refinement technique introduced by Colella and Woodward in |WC84] . 
We determine the interaction region is not only contained in the cone of sound, but it is 
exactly synchronized with the cone of sound region. We also show convergence of the non- 
interaction regions, the FRW-1 and TOV metrics, by computing the numerical error against 
their true solutions. We resolve the question of the nature of the secondary wave in this 
solution as being another shock wave reflected back in towards the center. We show the 
discontinuities in the fluid variables and metric derivatives, confirming these solutions to the 
Einstein equations are solutions in the weak sense of the theory of distributions. Changing 
the parameter f results in a quantitatively different solution; thus, we have a one parameter 
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family of quantitatively different shock wave solutions to the Einstein equations. The FRW- 
2/TOV matched model simulation gets setup in SectionEHwith the results shown in Section 
I7.5[ We numerically show the FRW-2/T0V matched model is the same solution as the 
FRW-l/TOV model, differing by a non-linear time coordinate transformation. This second 
matching provides a pedagogically interesting numerical confirmation of the covariance of 
the Einstein equations in standard Schwarzschild coordinates. 



7.1. Preliminaries 

In our simulation, we track the outer boundary of the lightlike and sound like information 
emanating from the initial discontinuity and compare them to the trajectory of the shock 
waves. In order to make this comparison, we need to determine how fast light and sound 
propagates through our universe. Since space and time are bent differently across the 
universe, we expect these speeds to be dependent on the point of interest. We first compute 
the (coordinate) speed of light. To simplify the calculation, we suppress the f'^dQ'^ term 
and take the 2-dimensional metric in standard Schwarzschild coordinates 

(7.2) ds^ = -Bdf + A-^dr'^ = gijdx'dx^. 

Consider a lightlike curve, 7(^) = (t(^),f(^)), which has the velocity vector 

Since 7(^) is a lightlike curve, the length of this curve is zero or 

(7.4) „.,,x*._«(|)^^-.(l)^ 

Solving this equation for the coordinate speed, we obtain the speed light travels as 

(7.5) t . I = ±VaB, 

taking the plus or minus sign for light traveling in the positive or negative direction, re- 
spectively. 

To find the speed at which sound travels in (f, f)-coordinates, we first note the sound 
speed is tied to the equation of state p = ap as being ^/a. This speed is relative to the 
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background fluid at rest (i.e. v = 0). To find the speed of sound relative to fiuid moving 
at a speed v, we need to apply the Lorentz transformation law for velocities (14.24p . More 
specifically, the speed ^/a transforms into 



(7.6) 



1 -L "V^ ' 



Again the plus/minus sign determines the direction in space being traveled. In order to 
determine the {i, f)-coordinated speed, we magnify this speed by the factor VAB, and the 
speed of sound becomes 

V ± y/a \ 



(7.7) s± = VAB 



1 ± "v°" 



Notice how the speed of sound is not only dependent on the metric {A,B), but also on 
the movement of the fluid v, indicating the dependence on the medium the sound travels 
through. 

For our simulation, we have a fixed speed of sound, namely ^/a = c^/\/3, and based 
on the position in space, we know the values for (^4, B) and v to determine both the speed 
of light (|7.5p and the speed of sound (17.7p . Hence, the new position of the light/sound 
information f is determined from the old position f^, after a time change of At by 

(7.8) r = n + sAi, 

where we set the speed as s = l± or s = s± for the light or sound information, respectively, 
choosing the plus/minus sign depending on the direction traveled. 

For the cone of sound verses the shock wave position comparison and error estimates, we 
want to know where the FRW metric stops and the interaction region begins; therefore, we 
need a mechanism to determine the border between the two. Also, we need to find a similar 
border between the TOV metric and the interaction region. Theoretically, the position of 
the shock wave can be used as the borders, but since the solution is numerically simulated 
with a first order method, these shock waves are smeared by numerical diffusion. This 
diffusion bleeds into the metrics, causing numerical error in testing this region against the 
known model. In order to handle these errors, we develop a criteria, based on studying the 
numerical solution, to determine where the numerical diffusion ends and the metric begins. 
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The border criteria we choose for the FRW side is where the spatial derivative of the fluid 
velocity first changes sign. In terms of notation, the border is the first position f such that 



This criteria is sufficient because the velocity is an increasing function of f for the FRW 
metric in standard Schwarzschild coordinates and only decreases when it hits the shock. 
To implement this criteria, we start the gridpoint under consideration at the minimum 
radius Xi = xi = fmim which is in the FRW metric because of the boundary condition, and 
increase the index i until the numerical derivative changes sign, giving us the point where 
the numerical diffusion first takes place. 

Next, we explore the border criteria for the interaction region and the TOV metric. 
Again, the fluid velocity is used as our indicator. Since the fluid velocity is theoretically 
zero, we would ideally flnd the point where it first becomes non-zero, but in the simulation 
it is only approximately zero, due to numerical error. Instead, the border criteria we choose 
for the TOV side is where the absolute value of the spatial derivative of the fluid velocity 
first becomes greater than 0.01. More specifically, the border is the first position f such 
that 



This criteria detects the first significant change in the velocity and ignores small changes 
occurring from numerical error. To implement this strategy, we start the gridpoint at the 
maximum radius Xi = Xn = Tmax, which is in the TOV metric, and decrease the index i 
until the numerical derivative satisfies (j7.10p . 

To run our simulation, we need boundary data from both the FRW and TOV metrics 
to maintain a consistent solution at the edges of the simulated region of the universe. 
Much to our surprise, the TOV metric boundary condition changes although the TOV 
metric itself is time independent. Even though the theory predicts there is a TOV metric 
matched continuously to the right side of the interaction region, the TOV metric allows for 
arbitrary changes of the time coordinate t, and we do not know which of these is matched 
in the solution except at the initial matching by solving for the correct Bq. More precisely. 



(7.9) 




(7.10) 



dv , , 

-(f) > 0.01. 
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consider an arbitrary time coordinate map ip : t ^ t applied to the TOV metric 



(7.11) 



B{r)dt^ + 



1 - 



2gM(f) 



1 



r 



) 



df2 + f2(iSl^ 



Since this metric is time independent, the effect of t = ^{t) is just a scahng of the B metric 
component, resulting in the transformed metric of 



In (r, f)-coordinates, the metric is no longer time independent because the time metric 
component now has a time dependent scale factor, but the other metric component along 
with the fluid variables are still independent of time. Since we refer to it throughout the 
chapter, we define 



to represent this scale factor caused by the arbitrary time coordinate transformation of the 
TOV metric. This allowance for arbitrary time transformations gives us a slew of potential 
TOV metrics, differing only in the time metric component, to match in our solution. One 
can view as a uniform change within the TOV metric of the rates the clocks move, 
and remarkably the FRW/TOV shock wave solution picks out the correct one of these TOV 
metrics at each time step based on the integration of B up through the FRW metric and the 
interaction region. To adjust to the change in the time scale, we determine where the TOV 
metric starts, based on the border criteria above, and match B^ at that border to determine 
which TOV metric is used as the boundary condition. This rematching also enables us to 
pick out the correct TOV metric model to test against the simulated solution in order to 
prove convergence of the TOV side of our universe, which will be discussed in more detail 
later. 

To explain this in different words, the theory in |JGB06] tells us the reason we must 
include the time scale factor in the B component of the metric is because in the numerical 
method the metric values on the left hand boundary of the simulated region are imposed, 
but the metric values on the right hand boundary are determined by the integration of 
the equations. Beyond the light cone on the right hand side, we know beforehand the 



(7.12) 




(7.13) 
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method has to simulate the TOV spacetime, but we do not know a priori in which standard 
Schwarzschild coordinate system it is simulated. Thus, since the standard Schwarzschild 
form allows for an arbitrary time scale factor on the metric B component, the correct time 
translation must be included to get the correct B component on the right hand boundary. 

7.2. FRW-1 and TOV Matched Model Setup 

In this section, we cover the details of building the initial profiles and the boundary data 
required to setup the FRW-1 /TOV model simulation. We restate both metrics for ease of 
reference. Using the coordinate transformation 

f =Vir, 

(7.14) r fM r2 



4^2] '4 

the FRW metric in standard Schwarzschild coordinates (the FRW-1 form) is 

(7.15) ds^ = -—^dP + —^df'^ + f^dn'^, 

1 — ■u^ 1 — 

with the metric components 

(7.16) A{0 = l-v{0\ B{0 



and the fluid variables 



(7.17) p(^,r) = -^, v{0 = ^ , 

where the variable ^ is defined as ^ = f/t and is a function of v by the following relation 

Whereas, the TOV metric is written 
(7.19) ds' = -B{f)dP + ( I df' + f'dn\ 



1 



with the metric components 



4ct- 



(7.20) A{f) = l- SttQj, B{f) = Bo{f) ^ , 
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and the fluid variables 



(7.21) p{f) = ^, vir) = 0, 



where the parameter 7 is a constant dependent on a 
(7.22) 7 ' ^ " 



27rg VI + 6(7 + 0-2 

Recall, for this simulation we let a = 1/3. 

Assuming an initial start time to > to be determined later, we choose an initial radius 
for the discontinuity of fo, and we want to match the metric components {A, B) continuously 
at the starting point (to,ro) for the initial discontinuity in the fluid variables. We start by 
matching the metric A component at this point 

(7.23) ApRwih, ro) = l-v (^S^ = 1 - St^Q-i = ATov{io,fo). 

Let vq = v{rQ/iQ) represent the fluid velocity at the discontinuity so (j7.23p implies 



(7.24) = 



4C7 



1 + 60- + a2 ' 

where we substituted (|7.22|) for 7. Take note that vq is independent of the our free parameter 
ro, it is quite astounding the velocity of the fluid at the discontinuity remains the same 
regardless of the placement of the discontinuity. Equipped with the value of vq, we use 
()7.18p . rewritten as 

to 1 + ^0 
to find the unknown starting time Iq as 

(7.26) to - • 

The independence of vq from fg along with (|7.26|) implies the initial start time is proportional 
to the initial radius of the discontinuity. Finding to enables us to build the initial profile 
of the FRW-1 metric for any radial coordinate f < fo by computing ^ = f/tg and using 
equations ffAI^ - fm) . 
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To compute the TOV metric, the A metric component is aheady determined beforehand 
by a constant (I7.20p . To find the other metric component, match it at the discontinuity 

1 



(7.27) 



forcing the constant Bq to take the form 
(7.28) Bo 



4o- 

" l + o- 



BFRw{to,ro) 



l-vl 



With the TOV time scale Bq, we can build the TOV metric for any radial coordinate f > fo 
by using the equations (j7.19p - (j7.2ip . 

Combining all this data together, we build the following functions Vinit{r), pinit{r), 
Ainiti^), and Binit{r) to use as the initial data at time Iq ()7.26p . depending on the free 
parameter tq. Because the other functions are based on v in the FRW-1 space, we start by 
stating the fluid velocity 







r < ro 
f > fo, 



(7.29) Vinitir) 
where ^ = r/tQ. Based on this function Vinit, the initial density profile is 
(7.30) 



Pinit{r) 



7 



r > fo. 



and the metric Ainu = (Ainu, Binit) is giving by 



(7.31) 



and 



(7.32) 



Hnit 



[r] 



Binit (r) 



1 — 87r^7 f > fo. 



r < ro 



l-vf. 



Bo{r) i+i^ r > ro. 



One last piece of information needed to run the simulation is the boundary conditions 
or the ghost cells. Recall, we need both boundary conditions to evolve according to the 
Einstein equations to maintain consistency at the edges of our simulated universe. For the 
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FRW-1 side, at the gridpoint xq, the fluid velocity is a function of time tj 

(7.33) voj = , 

where the variable ^ is defined as ^ = x^/ij. With the fluid velocity, the fluid density 
becomes 

(7.34) poj = 

Since the metric components are staggered relative to the fluid variables, we need to compute 
the half gridpoint 

Ax 

(7.35) xi=xq + — , 

2 2 

and use it to find the corresponding velocity 



(7 1- v^T^ 

(7.36) vy = , 

for ^ = xijtj. We use this velocity in the following computation for the metric components, 

(7.37) Ai,, = l-t-i B^, = -^. 

2'J I — Vi . 

The boundary condition for the TOV side is easier to implement. Since the TOV metric 
is independent of time, we set the data values during the initial setup of the function profiles, 
and they remain the same. These static values work for the fluid variables and the metric 
component A, but the function i?, specifically the time scale (j7.13p . changes and must be 
rematched during each time step, as discussed in the last section. Using the above criteria 
for the TOV border (j7.1Up . let Xi be the gridpoint position of this border. We rematch the 
time scale at time ij by the following formula 

(7.38) B'' = B{tj,Xi)[xi)-^, 

where B{tj,Xi) is the simulated solution at the coordinate (tj,Xi). 
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Figure 7.2. Initial profiles 



7.3. FRW-1 and TOY Matched Simulation Results 

In this section, we look at results of the simulation for the FRW-l/TOY matched model. 
We use the initial profiles (j7.29p - (|7.32p along with the boundary conditions (I7.33p - (l7.38p 
developed in the last section to run the simulation. Figure 17.21 shows the initial profiles 
with these parameters for the fluid variables {pinit^Vinit) and the metric P^mit-, along with 
the mass. By selecting the initial discontinuity at fg = 5, equation (I7.26P gives the initial 
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Figure 7.3. Evolution of the fluid variables during a unit of time 



start time of to = 5.4554. Notice how the discontinuities in the fluid variables jump down 
from the FRW-1 side to the TOV side. Moreover, the FRW density p and the TOV density 
p at this discontinuity are related by p = 3p. 

With these initial profiles, we run the simulation for one unit of time (i.e. tend = to + !)• 
Figure [7131 depicts the evolution of the fluid variables (p, v), giving us a frame by frame view 
for the evolution of the fluid variables across this time frame, evenly distributed from the 
left frame at to to the right frame at tend- After the initial time to, two shock waves are 
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formed, the stronger shock moving out toward the TOV side and the weaker shock moving 
in toward the FRW-1 side, creating a pocket of higher density expanding and interacting 
with both the FRW-1 and TOV metrics; therefore, the secondary wave to the strong shock 
for this solution of the Einstein equations is another shock wave, reflecting back in. 

Next, we focus our attention to the resulting solution at the end time, tend- Figure [73] 
highlights where the two shock positions are relative to the cone of sound and the cone of 
light. The cone of the light is represented by the white region while the cone of sound, 
embedded in the cone of light, is represented by the grey region. Notice how the edges of 
the cone of sound line up with both shock wave positions, showing the interaction region 
between the two metrics lie completely in the cone of sound. Since both characteristics and 
the edges to the cone of sound move at the speed of sound, we understand this result because 
these edges impinge on the shocks like a characteristic, so if one of the edges were to get 
slightly ahead or behind the shock position, then that edge would get pushed into the shock 
like all characteristics close to the shock. This figure also displays the spatial derivatives 
in the metric components A and B, the green and orange graphs, respectively. These 
derivatives {A',B'), found using numerical differentiation, have discontinuities aligned with 
the ones for the fluid variables at the edges to the cone of sound. Looking back at Figure 
17. 3^ it shows the profiles for the metric (A, B) as being continuous, so the metric is Lipschitz 
continuous, reinforcing the fact that we have a weak solution to the Einstein equations. 

The convergence of this solution, tested by successive mesh refinements, is shown in 
Table 17.11 Since we are implementing a first order method, we expect a convergence rate 
of 0.5 for the discontinuous fluid variables, while we expect a convergence rate of 1 for the 
continuous metric. Looking at the results in the table, the convergence rate of the fluid 
variables start around 0.5 for the big mesh sizes, as expected, and improve to 1 as we 
mesh refine, which is a little surprising. In contrast, the convergence rate for the metric 
components start around 1 and on average continue to stay around 1 under mesh refinement, 
regardless of the high/low swing in the convergence rate for the metric B component. We 
believe this change in the rates is due to numerical error from integrating the metrics up 
across the universe. 

Next we study the preservation of the FRW-1 and TOV metrics outside of the interaction 
region. Figure 17.51 shows the numerical solution after one unit of time for a mesh with less 
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Figure 7.4. Solution after a unit of time, showing the derivatives of the metric 



gridpoints {n = 1024). In this figure, the numerical diffusion is more pronounced as opposed 
to a finer mesh, as seen in Figure 17.41 This smearing of the shock wave is the reason for 
creating the border criteria discussed earlier. This figure also highlights why each border 
criteria works. The FRW border criteria is sufficient because one can see the fluid velocity 
on the FRW-1 side is increasing until it hits the smeared shock wave and starts decreasing. 
On the other side, the TOV border criteria is sufficient because the fluid velocity is almost 
constant until the smeared shock wave causes a significant change. Even though it looks 
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4.334e-004 
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2.568e-004 


0.76 


5.160e-004 


1.9 
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0.92 
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4.172e-004 
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1 
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1 
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16384 
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1 
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1 
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1.1 
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0.12 



Table 7.1. Successive mesh refinement convergence results 



constant on tlie TOV side, the fluid velocity is only close to being constant because there are 
numerical errors too small to be displayed on the graph shown. These borders are displayed 
in the simulation by where the grey lines stop, as shown in both Figure [731 and Figure EH 
Take notice on how both borders stop at the edge of the smearing of the shocks, and as 
we mesh refine, from Figure 17.51 to Figure 17.61 they tend towards the edge of the cone of 
sound, as desired. The colored lines verses the gray lines represent the simulated solution 
verses the model solution using equations (I7.15p - (l7.22p . We record the convergence between 
the FRW-1 border, based on the criteria above, and the left edge of the cone of sound in 
Table [721 Table [72] shows, from left to right, the number of gridpoints, the position of the 
left edge to the cone of sound, the position of the FRW-1 border, the error between the 
two, and the rate of convergence associated with the error. The corresponding results are 
displayed for the TOV side in Table 17.31 where we record the convergence between TOV 
border and the right edge of the cone of sound. Notice how the error between the cone of 
sound and both borders are decreasing, approaching zero. Moreover, the convergence rate 
is increasing, approaching 1, a linear rate which we would expect from a first order method. 
These results gives us confidence in our method for determining the edges to the cone of 
sound and the FRW and TOV borders. 

Equipped with a method for determining the position of the FRW-1 and TOV borders, 
these borders are used to show convergence of the metrics outside the interaction region. 
These borders are indicators on where to stop computing the error between the simulated 
and model solutions. For example, we consider the FRW-1 side error computation for 
n = 1024. Using Table [721 t^is numerical diffusion on the FRW-1 side ends at = 4.8495, 
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Figure 7.5. Smearing of the shock waves for solution with less gridpoints 
(n = 2^0 = 1024) 



so when computing the error, we only consider mesh points Xj such that xi < Xi < r^,. As 
we mesh refine, the point r* increases, getting closer to the cone of sound as shown in Table 
7.21 and the region where we are performing the error calculation expands, converging to 
the edge of the interaction region. Using this procedure, the results for the convergence for 
the FRW-1 and TOV side are listed in Tables [7^ and [731 respectively. These tables reveal 
both sides are converging to their respective model solution at the expected rate of 1. We 
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Figure 7.6. Showing the model metrics against the simulated solution 



find it remarkable that the TOV metric is preserved, regardless of all the mix-up in the 
interaction region. In particular, the TOV metric components remain intact even though 
we are integrating through the interaction region to obtain the metric. 

Since our free parameter is the position of the initial discontinuity, we have a one 
parameter family of solutions to the Einstein equations based on the initial radius, and 
we briefly explore changing this parameter. Figure 17.71 shows the initial profile for shock 
position at fg = 95, with fmin = 93 and fmax = 97. We make a couple of observations 
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Table 7.2. Shock wave verses cone of sound results for FRW side 



Gridpoints 


Cone of Sound 
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Rate 
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Table 7.3. Shock wave verses cone of sound results for TOV side 
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Table 7.4. Convergence results for the FRW side 



of the initial profiles, comparing to Figure 17.21 The fluid velocity along with the metric 
components at the initial discontinuity are the same, as determined in previous analysis. 
All these profiles are just stretched out over a longer region of space, causing them to look 
more like straight lines. There is less density and greater mass in this region of the universe 
since we are farther from the center of it. Figure 17.81 shows the solution after one unit of 
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Table 7.5. Convergence results for the TOV side 



time has passed, where again we have two shocks waves bounding a high density region. 
Notice for this solution the weak shock has a positive speed as opposed to the former case, 
as seen in Figure 17.41 The difference in speed from the ear her simulation fo = 5 indicates we 
are dealing with a quantitatively different solution from before. We explore changing this 
parameter many times and determine two conclusions in each case. One is the resulting 
solution always has a region of higher density surrounded by two shock waves, a strong 
shock on the TOV side and a weak shock on the FRW-1 side, and the other is the shock 
waves have different speeds, resulting in quantitatively different solutions. Hence, we truly 
have a one parameter family of quantitatively different shock wave solutions to the Einstein 
equations. 



7.4. FRW-2 and TOV Matched Model Setup 

After exploring the FRW-l/TOV model simulation, we cover the details in the setup of 
the FRW-2/T0V model. We recall the FRW-2 metric, so using the coordinate transforma- 
tion 

f =Vtr, 

(7.39) ^- 

the FRW metric in standard Schwarzschild coordinates (the FRW-2 form) is 

(7.40) ds"^ = —pTjT^ i^df + -^dr^ + f^dn^, 
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Figure 7.7. The initial profiles for the shock radius at ro = 95 



with the metric components 



(7.41) 



A{t,r) = l- v'^, B{t,r) 



1 



^2(1 _^2^ 



the integrating factor 



(7.42) 



t 

4t^ + f2 ' 
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Figure 7.8. Solution for the shock radius at vq = 95 after a unit of time 



and the fluid variables 

(7.43) p{i,r) 



3 — Tj 7* 

4^' = 2 = 2^- 



Remember, unlike FRW-1, the FRW-2 metric relies on the FRW time coordinate t, which 
is the following function of (t , f) 



(7.44) 



2^1 
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We match this metric to the TOV metric (equations (|7.19p - (j7.22p ). and we follow a proce- 
dure similar to the FRW-l/TOV matching. Again, we choose the initial shock position fo, 
assuming an initial start time to > 0, and we want to match the FRW-2 and TOV metrics 
continuously at the point (to,fQ). Like before, we match the metric A component as in 
equation (j7.23p to obtain 



(7.45) vo = ^/8^ ' 



1 60- a2 ^ 



where vq represents the initial velocity at the interface. Unlike before, we do not have a 
relationship between and to (|7.18p to solve for the initial time to- Instead, we possess 
a relationship between vq and the FRW time coordinate to in equation (j7.43p to find the 
FRW time coordinate at the discontinuous interface 

(7.46) to(to,fo) = ;^. 

2^0 

In order to find the initial time to , we solve for the integrating factor constant ^0 > requiring 
us to determine the integrating factor. In the FRW-1 matched model, the integrating factor 
constant was suppressed, which is equivalent to setting ^fo = 1, causing the (coordinate) 
speed of light \J AB on the FRW-1 side to be one. We follow the same paradigm with the 
FRW-2 matched model by choosing ^'o such that the speed of light is one on the FRW-2 
side, which is equivalent to setting = 1. More specifically, we choose \I'o such that at the 
discontinuity 



implying the integrating factor constant must be 



(7.48) *c 



/4t2 + f2 



V *o 

Now we use (|7.39p to solve for the initial time with fo, to (!7.46p . and (|7.48p as 



One might worry ^ does not equal one across the FRW-2 region for other values of f besides 
fo- This concern is not an issue since we can substitute the time coordinate ()7.49p into the 



7.4. FRW-2 AND TOY MATCHED MODEL SETUP 



110 



integrating factor equation (j7.42p to obtain 

(7.50) = = 1, 

Zto 

and to is independent of the spatial variable f. Equipped with to and ^'o, we are capable of 
solving for t(to,f) for any f < fo enabling us to compute the fluid variables and the metric 
components for the FRW-2 region. 

For the TOV metric, the A metric component is the same constant as the FRW-1 case 
()7.20p . so we are left with matching the B metric component. Matching BpRw and Btov 
at the coordinate (to,^o) provides us with 

- _ _ 4n- 1 _ _ 

-^o) 

Since ^' = 1, the time scale constant Br, becomes 



4o- 

■ 1 + CT 



2 ' 



(7.52) Bo ^ 

exactly the same as the FRW-1 case (j7.28p . With Bq, we can compute the TOV metric for 
any radial coordinate f > fo by using the equations ()7.19p - ()7.2ip . 

With the matching complete, we build the functions Vinit{r), Pinit{r), Ainit{r), and 
Binitir) to use as the initial data at time to. These initial profiles are similar to their 
FRW-1 matching counterparts (|7.29p - (|7.32p . except the velocity is a function of rj instead 
of ^, and we rely on the FRW time coordinate t{t,f) defined in (j7.44p . Again, due to the 
dependence on the function v, we state the fluid velocity first 



(7.53) Vinit{r) = 
where t] = f/t. The initial density profile is 

(7.54) Pinitif) = < 



^ r < ro 
f > fo. 



^ r > fo. 
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Based the function Vinit, the metric components are giving by 
(7.55) Ainitir) = 



1 - vfnit r<ro 
1 — SttG^ r > fo, 



and 



(7.56) Binitir) 



init > 
4a 



Bo{r) i+<^ r > ro. 



Since the boundary conditions are derived from the model equations, they contain the 
same discrepancies between the FRW-1 and FRW-2 matchings as the initial profiles. For 
the FRW-2 side, at the gridpoint xq at time tj, we use ()7.44p to obtain the FRW time 
tj{tj,xo) for computing the fluid velocity and density. We find the fluid velocity 

(7.57) vo,j = 
and the fluid density 

(7.58) poj 



Since the metric components are staggered relative to the fluid variables, we need to compute 
the half gridpoint as before ()7.35|) and use it to find the corresponding velocity 

(7.59) v.^^ = I 

for r] = xi/tj. We use this velocity in the following computation for the metric components, 

1 



(7.60) ^i. = l-l,' ^^^■ = ^2(1^^ 



Since the TOV metric is independent of the FRW metric matched with it, the TOV 
boundary condition remains the same, and we use the same matching as in the FRW-1 case 
(I738D . 

7.5. FRW-2 and TOV Matched Simulation Results 

Using the initial profiles and boundary conditions developed in the last section, we run 
the simulation of the FRW-2/T0V matched model, comparing it with the FRW-l/TOV 
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Figure 7.9. Initial profiles 



model results. Figure [7^ shows the initial profiles for all the variables in the FRW-2/T0V 
matched model. Comparing with the FRW-1 case (Figure [7.21) . the graphs in both figures 
match exactly. Actually, the only difference between the two models is the start time, 
to = 10.9109 instead of the previous time of to = 5.4554. This similarity is one indication 
that we might have the same solution as the FRW-1 case. 
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After running the simulation for one unit of time as shown in Figure I7.10| we notice 
more similarities to the FRW-l/TOV model. More precisely, if we compare with the FRW- 
1 case in Figure 17.61 both solutions have very similar features. These features include the 
formation of two shocks with the stronger shock moving outward and all five graphs share 
the same shapes. The main difference between the two solutions is the FRW-2 solution 
appears to have run for a longer period of time. This difference is highlighted by the light 
and sound like information traveling slightly farther out from the initial discontinuity. For 
example, the ingoing edge to the cone of sound is 4.8695 and 4.8576 in the FRW-1 and 
FRW-2 solutions, respectively, meaning the weak shock wave in the FRW-2 solution moved 
slightly farther out than the corresponding wave in the FRW-1 solution. This outcome 
suggests we have the same solution but at different stages of their development. Since both 
these solutions share the same initial data and the same TOV boundary condition, the 
only difference is the FRW boundary condition. Looking deeper and comparing equations 
(|7.33p - (|7.37p to ()7.57p - (|7.60p . if we assume the same FRW time t at this boundary, the 
only change in the boundary condition is the metric component B due to the presence of 
a non-constant integrating factor. Intuitively, it makes sense to have the same solution at 
different times because a change in B affects the measuring of time verses space. This does 
not affect the interactions just how time and space are bent relative to each other, and 
since we fix spacial distances by setting constant grid cell sizes, it affects how time changes 
relative to each grid cell. Hence, there is more evidence to support the claim both solutions 
are the same, just the clocks are not synchronized to produce the same result. 

We proceed under this assumption, a same solution at different times, and search for the 
FRW-2 time t2 where the solution at this time matches the solution of the FRW-l/TOV 
model at ti = 6.4554, running for one unit of time. If our assumption is correct, both 
models are mapping over the same region of spacetime in the original FRW metric. For 
this mapping to work, both times ti and t2 would have to correspond to the same FRW 
coordinate time t. Of course, since t is a function of f as well, this correspondence must 
hold for every f in the FRW region of the universe even though each f corresponds to a 
different time t. To find i2, we pick a radius to match the time t . We use the left boundary 
radius f^m and find the corresponding FRW-2 time to be t2 = 11.8688. This time is 0.9579 
units of time after the initial start time of to = 10.9109, as observed in Figure [7^ and this 
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Figure 7.10. Solution after a unit of time 



change in time is barely less than the change of one unit, explaining our earlier observation 
that the FRW-2/T0V solution is at a slightly later stage of development after one unit of 
time. 

After finding this corresponding time ^2, we run the FRW-2 model for 0.9579 units of 
time and show the results in Figure [7.111 Comparing against Figure [7l6l the solutions are 
almost exact, with the same shape and similar numbers, except for the B metric component. 
Interestingly enough, the metric B component graph has the same shape but different 
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Figure 7.11. Solution after 0.9579 units of time 

numerical values, so if both graphs are placed on top of one another, they overlap, meaning 
the metric component in the FRW-1 case is stretched and shifted relative to the FRW-2 
case. 

We continue on and test our same solution hypothesis numerically. We accomplish 
this task by using the FRW-l/TOV solution at a fine mesh refinement, n = 16,384, as 
the model to compare against. We perform convergence calculations assuming the FRW-1 
matched solution is the true solution, so we run the FRW-2/T0V matched model for 0.9579 
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Table 7.6. FRW-2/T0V model verses FRW-l/TOV model convergence results 



units of time for different mesh sizes and compute the error between this solution and the 
FRW-l/TOV solution, testing the convergence rates. One problem with this approach is 
the B metric components between the two solutions are shifted and stretched from one 
another, so a direct error calculation for B does not result in convergence. To overcome 
this discrepancy, we build a map between the two metric components, which we label as 
Bi and B2 for the FRW-l/TOV and FRW-2/T0V models, respectively. This map is built 
using the following data: let Bj,^^^ and B^^^ represent the maximum and minimum values, 
respectively, for the B metric component in the FRW-l/TOV model, and let -B^az ^^'^ 
be the corresponding values for the FRW-2/T0V model. We define the mapping 
from Bi to B2 as follows 

(7.61) B2 = f ^r^I^r ) {B^ - Bi,J + Bl,^. 

\ ^max min / 

The value -B^m represents the value the left most point gets mapped to while the scale 
factor (S^ax ~ Bmin)/(.-Brnax ~ ^min) represents the stretching. For the mapping between 
the solutions in Figure ESI and Figure Eni is 1.2598 while the scale factor is 1.1833. 

Using this map, the convergence results are shown in Table 17.61 ^-nd the data shows the 
FRW-2/T0V model is converging to the FRW-l/TOV model at a linear rate, confirming 
our same solution at different times hypothesis numerically. 

We find the discovery that the two FRW transformations produce the same result too 
uncanny to be just a random coincidence, and we search for a theoretical connection to 
justify this discovery. Since we are dealing with a solution to the Einstein equations of a 
spherically symmetric metric in standard Schwarzschild coordinates, we take another look 
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at the general form for this metric 

(7.62) ds^ = -Bit, f)dP + —L-df^ + r^dn^. 

^ ^ ^ ^ A(t,f) 

Let us consider the coordinate freedoms that preserves this metric; it turns out the only 
freedom of this metric is an arbitrary nonlinear transformation of the time coordinate t. We 
cannot change the radial coordinate f because the spheres of symmetry will no longer match, 
and we cannot change the time coordinate i in terms of f because this type of transformation 
would yield mixed terms violating the metric being in standard Schwarzschild coordinates. 
Therefore, if the spacetime is in standard Schwarzschild coordinates, the only change of 
coordinates of this metric, using the same center, is an arbitrary nonlinear change in the 
time coordinate. Since all the mappings from FRW into standard Schwarzschild coordinates 
differ by how the time coordinate t(t, f) gets mapped, based on different integrating factors, 
all these mapped FRW metrics must be the same solution which differ by a nonlinear 
change in the time coordinate. Hence, both metrics FRW-1 and FRW-2 developed earlier 
are actually the same solution with a nonlinear change of time between the two, implying 
the FRW-l/TOV and FRW-2/T0V matched models are the same too. 

With this fact in mind, we look for the transformation between the two time coordinates 
ii and t2 for the FRW-l/TOV and FRW-2/T0V models, respectively. Considering an 
arbitrary coordinate {t,f), the FRW-1 time ti is 

and the FRW-2 time t2 is 



(7.64) f, = i5, 

where we set the integrating factor constant 



(7.65) ^o = ^l!!±l!, 

as before (j7.48p to make = 1 and initially match the B metric component for the FRW- 
l/TOV and FRW-2/T0V models. Studying these equations, we notice a couple of things. 
The first one is the integrating factor constant chosen can be written as an explicit function 
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of h 



(7.66) 



^0 



4*2 + f2 



= 2 



4*2 + f2 




t 



Not only that, but the FRW-2 time can also be written as an explicit function of ti 



Looking at the specific coordinate {tQ,ro) for the initial discontinuity, we combine the equa- 
tions (j7.66p and (|7.67p to determine 



explaining the relationship between the FRW-l/TOV start time to = 5.4554 and the FRW- 
2/TOV start time to = 10.9109. We also build the relationship between the two end 
times using (j7.48p as the integrating factor constant in the time mapping equation (|7.67p . 
Substituting the FRW-l/TOV end time h = 6.4554 into this equation ^L67\\ gives us the 
corresponding FRW-2/T0V end time as t2 = 11.8688, matching the end time above to 
produce Figure 17.111 and to obtain the convergence results in Table 17.61 

The other mystery on our hands is the shifting and stretching of the metric B com- 
ponent. The shifting is a result of the difference in the boundary data of B between the 
FRW-1 and FRW-2 cases, caused by the different integrating factors. This metric compo- 
nent is computed at the end of each time step by integration of the solution up from the 
boundary, and a shift in the boundary data of B will result in a shift of the entire function 
B. The stretching is caused by the effect on the metric of the coordinate transformation 
from the ti to the t2 coordinate system. More specifically, consider a time coordinate map 
: t2 — )• ti from the FRW-2/T0V to the FRW-l/TOV time coordinate. This map produces 
the following relationship between the differential forms 



(7.67) 




(7.68) 



t2 = 2ti, 



(7.69) 
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causing a scaling or stretching between the B metric components. For our particular case, 
the time coordinate map is 

(7.70) h = ^{h) = , 

giving us a scale factor of 



(7.71) 




Plugging in the values for ^'q and the scale factor equals 1.1833, matching the observed 
value derived numerically above. This analysis gives us a procedure of mapping the FRW- 
1/TOV model over to the FRW-2/T0V model for any FRW-1 time ti, confirming our 
hypothesis that both models represent the same solution at different time coordinates. 
Consequently, running the FRW-2/T0V simulation is just an exercise in the time coordinate 
invariance of this solution of the Einstein equations in standard Schwarzschild coordinates. 



7.6. Summary 

We provide the specific details in building the initial profiles (shown in Figure 17. 2p and 
the boundary data required to run the simulation using the locally inertial Godunov method. 
We have the first glimpse of a shock wave in general relativity by providing snapshots of the 
solution in Figure [731 This solution results in two shock waves, a strong outgoing wave and 
a weak incoming wave, enveloping a region of higher density. Changing our one parameter 
f produces quantitatively different solutions, shown in Figure 17.81 This result enables us 
to resolve the secondary wave in the solution as a reflected shock wave. This shock wave 
solution is a weak solution to the Einstein equations because of the discontinuities in the 
derivative of the metric shown in Figure 17. 4[ where these discontinuities are aligned with 
the corresponding discontinuities in the fluid variables. We show convergence of the region 
of interaction to the cone of sound by recording the convergence between the position of the 
shock waves and the edges of the cone of sound in Tables 17.21 and 17.31 This result affirms 
that the region of interaction is precisely the region of the cone of sound, and again at this 
stage we have no formal mathematical proof of this. The numerical determination of these 
borders help us demonstrate the numerical convergence of the FRW and TOV metrics in 



7.6. SUMMARY 



120 



the non-interaction regions, using the exact solutions, in Tables [71^ and [731 because these 
borders gives us the boundary of integration in the calculation of the one-norm. 

In simulating the FRW-2/T0V model, we force the initial profiles for this model to 
closely resemble the FRW-l/TOV model, and after running the FRW-2/T0V model, we 
notice many similarities to its counterpart, providing evidence that these are the same 
solution at different times. We test and confirm this hypothesis numerically by proving 
convergence between the two solutions in Table [7l6l We proceed to find the theoretical jus- 
tification of these numerical results and determine the two solutions differ by a non-linear 
time coordinate transformation. This result can be interpreted as a numerical confirmation 
of the covariance of solutions to the Einstein equations in standard Schwarzschild coordi- 
nates. All these results give us confidence in the correctness of our locally inertial Godunov 
method implementation and in the accuracy of our solution. 



121 



CHAPTER 8 

Time Reversal Model 

Now that we have explored the shock wave model by simulating the FRW/TOV matched 
model forward in time, we consider reversing time and running the matched spacetime 
backward in time. Not only are we interested in the structure of the resulting solution, but 
we conjecture reversing time will send the strong wave into the center of the universe in 
hopes of forming a black hole. By sending matter into the center, we expect the density 
and consequently the mass function to increase as time unfolds. This increase at a fixed 
radius f will cause the black hole criteria 

(8.1) 2eMfir)^^_ 

r 

to eventually be satisfied for some time-space coordinate (t, f) . Since it is referred to 
throughout the chapter, we define the black hole number fj,(t, r) as 

/ X X 2gM(i,r) 

r 

Section 18.11 discusses the validity of a reversed time solution along with the setup to 
implement it. We show the difference of this setup from the forward time solution is just 
the sign of the times, i and t, from positive to negative, so the initial profiles and boundary 
conditions developed in Section 17.21 are used again here with this slight change. In Section 
18. 2^ we run this reversed time simulation for one unit of time to obtain results of the 
solution. In particular, we demonstrate numerical convergence of the entire solution, the 
interaction and non-interaction regions, and within the non-interacting regions, numerical 
convergence to the pure FRW and TOV metrics is also verified. As in the forward time case, 
we demonstrate the region of interaction is synchronized with the cone of sound region. This 
solution contains two expansion waves surrounding a region of under density, as opposed 
to the two shock waves in the forward time model. Both of these waves are an expansion 
wave in curved spacetime which we interpret as a generalized rarefaction wave, and we refer 
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to it throughout as a rarefaction wave. Since we are dealing with rarefaction waves, the 
metric is continuous, producing a strong solution to the Einstein equations. After one unit 
of time, the strong rarefaction wave is heading toward the center of the universe, bringing 
mass along for the ride, and getting closer to the black hole criteria (|8.ip . In Section \8.3\ 
we continue the flow of time to see the beginning of what we believe is the formation of a 
black hole, where we obtain a black hole number of 0.9218. Although we wish to get even 
closer to satisfying the black hole criteria (j8.ip . we are content with passing the Buchdahl 
stability limit of 0.9, which loosely states, after reaching this limit, a star is unstable and 
is subject to it own gravitational collapse. That is, the theorem of Buchdahl is whenever 
a mass gets within 9/8ths of its Schwzarzschild radius, or in our language its black hole 
number gets within 0.9, there doesn't exist a static solution of the TOV equations with 
a finite pressure at the center capable of holding the star up |ST97bj . This fact gives us 
confidence that not only does black hole formation occur in the reverse FRW/TOV model, 
but our simulation gives us a first glimpse at it. 

8.1. Reversing Time 

We argue the time reversibility of the FRW metric. Recall, the FRW metric is given by 

(8.3) ds'^ = -dt'^ + R^{t)\-—^—.dr^ + r^dn 



Plugging (18. 3p into the Einstien field equations (12. ip . along with the perfect and comoving 
fluid assumption, gives the following pair of constraint equations on the functions R{t), p{t), 
and p{t): 

(8.4, 

(8.5) R^ + k= ^(,R\ 

o 

developed by Smoller and Temple in jST95]. For the case of interest, when the metric 
is spatially flat [k = 0) with an isothermal equation of state {p = ap), these constraint 
equations ()8.4p and ()8.5p simplify to 

(8.6) '^ = T'^' 
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(8.7) = ^-fpR\ 

There exists explicit solutions, p{t) and R{t)^ to these pair of ODEs for t > |ST95) . 
Knowing these solutions exist, we show p{—t) and R{—t) for t < 0, the reversed time 
solutions, are also solutions to these ODEs. Define r = — t, and since r > there exists 
solutions p(t) and R{t) to ()8.6p and (j8.7p . Using these solutions and the fact ^ = —1, we 
show 

. dtdp (CT + I)p p (a + l)pdi?dt (cT+l)P p 

0.8 — P = -, — T = Pt = — Rt = ;^ = ;^ — -K, 

^ ' ^ drdt ^ m 3R dt dr 3R 

implying 

(^•^^ P - 3R{-t) 

The other ODE is easily satisfied by 

(8.10) R^ = i^^^kj = Rl = '-fpR^ = '-fpi-t)R\-t). 

Hence, the time reversed solutions p{—t) and R{—t) for t < are also solutions for the FRW 

metric. Notice the FRW metric is invariant under this time reversal r = —t: 

(8.11) 

ds"^ = -dT^ + R^{t) | ^_^^^^ dr^ + r^dJ^^I = -df + R^i-t) | ^_^^^^ dr^ + r^dn"^ 

We turn our attention to how this time reversal affects the FRW-1 metric. If we re- 
peat the process in Chapter [3] to obtain the coordinate transformation, the time reversed 
transformation becomes 

(8.12) ( -2 ^ 



4t\ 

and the resulting equations for the fluid variables (p, v) and the metric components {A, B) 
remain the same. The only difference is that t < implies t < and consequently, ^ = f/i 
changes sign from positive to negative. This change in sign only affects the fluid velocity 
since all the other equations in the metric posses f ^ terms in them. With this in mind, we 
use the results in the FRW-l/TOV matched model process in Chapter [7] with a negative 
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Figure 8.1. Initial profiles 



time t < 0. In particular, we use the same equations (I7.29p - (l7.38p for the initial profiles 
and the boundary conditions with a negative time, increasing as the simulation runs. The 
only difference is the fluid velocity at the discontinuity vq is chosen to be negative instead 
of positive. 
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Figure 8.2. Evolution of the fluid variables during a unit of time 



8.2. Reverse FRW and TOV Matched Simulation Results 

To better compare and contrast the forward and reverse solutions, we use the same initial 
parameters as before (17. ip . Figure [8T] shows the initial profiles for the fluid variables, the 
metric components, and the mass function. Comparing these profiles to the corresponding 
ones for the forward time in Figure 17.21 we see they match except the fluid velocity along 
with the start time = —5.4554 has changed signs, as expected from the above analysis. 
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After setting up the initial profiles, we run the simulation for one unit of time, but since 
time is running in reverse, to is getting closer to zero, evolving toward the big bang instead 
of away from it (i.e. tend = ^ o + 1 = —4.4554). Figure [52] gives a frame by frame view 
for the evolution of the fluid variables across this time frame, evenly distributed from the 
left frame at tQ to the right frame at tend- From the initial discontinuity, two rarefaction 
waves are formed, a stronger wave moving in towards the FRW region and a weaker one 
moving out towards the TOV region, and between the two waves a pocket of lower density 
is formed. These results differ from the forward time case where there are two shock waves 
surrounding a region of higher density. As time progresses, the rarefaction waves and the 
pocket are both expanding. To the left of this pocket, there is a density spike growing and 
moving in towards the FRW region, meaning more matter is falling into the center of the 
universe as time progresses. One can notice that indeed the FRW region is accumulating 
more mass from the start time (Figure l8.ip to the end time (Figure l8.2p . If more mass 
keeps coming into a finite radius, then fj, will continue to grow and eventually the black hole 
criteria will be satisfied ()8.ip . leading us to believe, given enough time, a black hole will 
evolve out of this solution. We explore this idea in more detail in the next section. 

We turn our attention to the end result at time t^nd, where Figure [8^ gives us a more 
detailed look at the fluid variables. This figure allows us to see the rarefaction waves in 
their entirety, where one might mistake them as shock waves in Figure 18.21 Comparing to 
the forward time model (Figure [7. 4p . the cone of sound is wider and shifted more to the left, 
while the cone of light remains the same. This shift toward the FRW side is caused by the 
negative velocity forcing the fluid and consequently all sound like information inward instead 
of outward. Observe the derivatives of the metric components are continuous, producing 
no jump discontinuities as seen in the forward time case; without these discontinuities, we 
no longer have to settle for a weak solution, so the reversed time model is a strong solution 
to the Einstein equations. 

Again, we use the successive mesh refinement technique to test for the convergence of 
this solution, recorded in Table [HTTl Since we are dealing with a continuous solution using a 
first order method, we expect convergence rates around 1 across all the variables. Looking 
at Table 18.11 the only variable close to this rate is the metric component A, which has an 
average convergence rate of about 0.9. The fluid variables have a rate close to 0.5, while 
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Figure 8.3. Solution after a unit of time, showing the derivatives of the metric 



the other metric component has a rate close to 0.7. It is unclear why these rates are lower 
than the forward model (Table 17. ip , where the discontinuities of the shock waves should 
make these rates lower, but numerical convergence is obtained at an adequate rate. 

We examine the preservation of the FRW and TOV metrics outside the region of in- 
teraction. Looking back at Figure 18.31 the ends of the rarefaction waves are outside the 
cone of sound. Before, there was numerical diffusion of shock waves bleeding outside the 
cone of sound, but we do not expect to see a similar phenomenon with a rarefaction wave. 
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Table 8.1. Successive mesh refinement convergence results 



Gridpoints 


Cone of Sound 


Shock Wave 


Error 


Rate 
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N/A 
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0.3 
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0.24 
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Table 8.2. Shock wave verses cone of sound results for FRW side 



Gridpoints 
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Table 8.3. Shock wave verses cone of sound results for TOV side 



Since the speed at which matter travels is bound by the speed of sound, this fact is a cause 
for concern. Like the forward model, we need to determine where the wave ends, and the 
border criteria developed in the Chapter [7] is used here. On the FRW side, the border 
criteria is where the derivative of the fluid velocity changes sign, and instead of increasing 
then decreasing in the forward time model, the fluid velocity is decreasing until it hits the 
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Figure 8.4. Showing the model metrics against the simulated solution 



rarefaction wave and increases, satisfying this criteria. On the TOV side, the border criteria 
is a significant change in the velocity from zero which is stih satisfied since the fluid velocity 
has a sudden decrease at the rarefaction wave, as opposed to the sharp increase at the shock 
wave in the forward time model. Looking at Figure ISTSl the borders are indicated by where 
the gray curves stop and shown to be at the end of both rarefaction waves, as desired. Next, 
we numerically compare the edges of both rarefaction waves against the cone of sound. The 
results for the FRW side are recorded in Table [8?2l while the TOV side is recorded in Table 
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Table 8.4. Convergence results for the FRW side 
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Table 8.5. Convergence results for the TOV side 



18.31 The convergence rates are significantly lower from before, but numerical convergence 
is still obtained, albeit a little slower. It is uncertain why both convergence rates are much 
slower, but these results give us confidence, as we mesh refine, that the outer edges to the 
rarefaction waves will converge to the edges of the cone of sound. 

Using these border criteria, we test the numerical convergence of the FRW and TOV 
metrics outside the region of interaction. Recall, these borders act as a marker where we 
stop computing the error between the simulated and model metrics, which is discussed more 
thoroughly in Chapter [71 Using this procedure, the convergence results for the FRW and 
TOV region are recorded in Tables [83] and ESj respectively. Looking at both tables, both 
regions are converging to their respective models at the appropriate rate of 1 for a first order 
method. Again, we note how remarkable it is the metric components in the TOV region 
remain preserved after integration through the FRW metric and the region of interaction. 
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Figure 8.5. Approaching the black hole with the interaction region hitting 
the TOV boundary 



8.3. Black Hole Formation 

In the previous section, as time progresses in the reversed time model more mass falls in 
towards the center of the universe. If this trend continues, then it is inevitable, given enough 
time, a black hole will form. In order to run the simulation long enough for a black hole 
to form, our initial parameters must be changed. Using the previous parameters ()7.ip . the 
region of interaction hits the boundaries before the black hole forms, and once this happens. 
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Figure 8.6. Looking at the black hole criteria and the speed of light 



the solution is no longer valid because our boundary data is inaccurate. In order to solve 
this problem, we expand the region of space under consideration by setting the minimum 
radius fmin = 0.1 and the maximum radius fmax = 20. We leave the other parameters the 
same, namely, the discontinuity position f o = 5 and the number of gridpoints n = 16, 384. 

Running this simulation, the interaction region hits the TOV boundary before the black 
hole forms, with the end result shown in Figure 18.51 We cannot proceed further because 
the TOV boundary data is no longer valid. At this stage of the solution, the highest value 
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Figure 8.7. Continuing black hole formation until it hits the boundary data 



of the black hole number is /x = 0.9110, located where the vertical magenta line is placed 
at a radius of f = 2.38. This radius corresponds to the dip in the velocity and metric 
component A. It also occurs right after the density spike which causes the mass to increase 
dramatically in this region of the universe. Comparing to an earlier stage of the solution 
in Figure 18. 4| at the edge of the incoming rarefaction wave, the density has continued to 
increase while velocity has become more negative, so matter has been coming in towards 
the center at a faster rate. Another thing to notice is the metric B component has grown 
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extremely large at the TOV edge of our simulation. This large metric B component causes 
a huge time dilation between both edges of our simulation, which can be seen more clearly 
in Figure where we graph the (coordinate) speed of light V AB in yellow and the black 
hole number in cyan. The black hole criteria is the most satisfied in a region of space 
around the point f = 2.38 where /u = 0.9110. Now more extreme relativistic effects come 
into play because the significant time dilation, measured by the speed of light, between the 
two edges is around 11,267 to 1, which means light travels 11,267 times faster at one edge of 
the simulated region to the other. This is a problem because it causes our future time steps 
to be smaller, so the evolution of the system with greatest fi is slower than the development 
close to the TOV border. As a final note, there is a blip in the velocity profile as seen in 
Figure 18.51 Since there is much data to support the accuracy of our numerical simulation, 
we have no reason to believe this blip is due to numerical error. We interpret this blip in 
the velocity as an unknown phenomenon, and it is a future research project to determine 
the nature of it. 

In order to get closer to the formation of the black hole, we could try to expand the 
region of simulated space again, but the huge time dilation causes the region of interaction 
to grow extremely fast, quickly reaching the new border while only a small amount of time 
passes near the region where the black hole appears to be forming (i.e. fi = 0.9110 at 
f = 2.38). Regardless of how far we put the maximum radius, we never get closer than 
/i = 0.9110 in the black hole criteria (j8.ip . Our alternative strategy is to realize we are only 
interested in the region of space where the black hole is forming, where fi is greatest. With 
this in mind, to continue the evolution of the black hole, we start chopping off the right side 
of our region of simulated space; even though the right most grid cell is no longer valid, the 
one next to it (on the left side) is still a valid solution. This process is like zooming into 
the region of the black hole formation. More precisely, the grid point Xn+i is used as the 
boundary data throughout the entire simulation. After the TOV edge has been reached, 
we use the gridpoint x„ as the new boundary data for the next time step, discarding the 
data associated with the gridpoint Xn+i- During the subsequent time steps, we discard the 
current boundary data at the gridpoint Xi and replace it with its predecessor Using 
this procedure and stopping when the radius of the right boundary matches the radius with 
the greatest value of the black hole number, we produce Figure 18.71 and Figure 18.81 The 
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Figure 8.8. Black hole criteria and the speed of light after the boundary 
data hits the black hole formation 



black hole number is /i = 0.9218 at a radius of f = 4.5, so we are a little closer to black hole 
formation. More interestingly, this radius for the greatest /i is moving farther out, from 
f = 2.38 to f = 4.5, meaning the region of black hole formation is expanding which is shown 
in Figure [8^ In Figure \87l\ one notices the metric component B, along with the speed of 
light in Figure ESI is growing exponentially at this radius, causing the time dilation to be 
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more extreme. Here, the time dilation from the center to this radius is 29,517 to 1, two and 
a half times greater than before. 

8.4. Summary 

We began the chapter by arguing for the validity of the existence of a reversed time 
solution, and then setup the initial profiles and boundary conditions needed to start the 
simulation. It turns out the only difference from the forward model is that positive time is 
replaced by negative time, and this difference only affects a change in sign of the velocity, 
as shown in Figure 18.11 Snapshots of the solution are recorded in Figure 18.21 showing the 
reverse time solution containing rarefaction waves surrounding a region of under density. 
Unlike the forward time model, the derivative of the metric is continuous in Figure 18.31 
so this model is a strong solution to the Einstein equations. We mimic the convergence 
results obtained in the forward time model. In particular, in Table \87i\ we record numerical 
convergence of the entire solution taken as a whole, including the interaction and non- 
interaction regions, using successive mesh refinement. In Tables (8^ and [8131 we also record 
numerical convergence of the interaction region to the cone of sound, using border tracking 
of the edges to the simulated rarefaction waves. In Tables 18.41 and 18.51 we record the 
numerical convergence of the non-interaction regions, using the true solutions. Like the 
forward time model, these results affirm the region of interaction is precisely the region of 
the cone of sound (no formal mathematical proof of this is known). 

As time progresses in this solution, more and more mass is falling in towards the center 
of the universe. This fact leads us to believe that if the solution were to progress in time 
beyond the simulation time limit, the black hole number (j8.2p would continue to increase 
and eventually the black hole criteria (|8.ip would be achieved exactly, resulting in black hole 
formation. Our method is to expand the region of simulation by running the simulation 
until the TOV region is completely gone, giving us a value of /x = 0.9110. We then continue 
to run the simulation by discarding the incorrect boundary data and zooming into the region 
of black hole formation, and by this refined method we obtain convergence all the way up 
to the black hole criteria value n = 0.9218. These results taken together give us confidence 
in our solution as a whole and provide strong support for our claim that this solution 
has black hole formation contained within it. As a final note, with the successful results 
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obtained here and Chapter [71 we have simulated general relativistic analogs of both the 1 
and 2 family of shock and rarefaction waves, the elementary waves of conservation laws. 
This fact gives us confidence in the ability and accuracy of our locally inertial Godunov 
method to simulate any solution to the Einstein equations for a perfect fluid in standard 
Schwarzschild coordinates. 
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CHAPTER 9 



Putting Units into the Simulation 



We end our discussion of shock wave formation by putting back the units to give more 
physical meaning to our simulations. To place units on our numerical values requires us 
to understand the consequence of setting Newton's gravitational constant Q and Einstein's 
speed of light constant c both to one. The effect of this is the unification of the units of 
time, length, and mass into one set of units, which we choose as the units of mass in our 
discussion here. Understanding this unification allows us to recover the other units, time 
and length. To end this chapter, we consider our simulation on the solar and galactic scale 
by transforming our numerical values to familiar units, giving us a firm grasp of the scale 
and magnitude involved within our simulations. 

To determine the effect of setting the speed of light constant c to one, we start by 
considering the TOV metric 



where we dropped the f^dQ^ to simplify the notation. Since it represents a measurement, the 
metric must have units of length or time attached to it; these units can be either associated 
with the coordinates {t, f) or the metric components themselves. Here, we choose the units 
to be associated to the coordinates (f, f), and we use units of length (i.e. kilometers) to be 
the measurement of the metric. This implies the time coordinate t is measured in units of 
length, and the speed of light constant converts units of time (i.e. seconds) into units of 
length. More precisely, let be the time measured in units of time, so the time coordinate 
becomes 



(9.1) 




(9.2) 



t = 



ct: 



where the units are represented by units of length scale. One can view this as the time 
coordinate absorbs the speed of light constant and transforms it into units of length. This 
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absorbtion is equivalent to setting c = 1, as done in our simulations. In particular, the 
following holds 

df df df 

(9.3 = = = 1, 

dt^ d[ct^) dt 

so the coordinate speed of light is scaled down to one in (t, f) coordinates. To recover the 
units of time from the time coordinate i requires a simple calculation 

(9.4) U = -. 

c 

Thus, the speed of light constant enables us to make equivalent the units of time and length, 
a truly remarkable insight of Einstein's theory. 

To determine the effect of setting Newton's gravitational constant Q to one, we perform 
dimensional analysis on (19. ID to determine the dimensions of ^. To this end, let L represent 
units of length, T represent units of time, and M represent units of mass. As stated before, 
the metric components of (j9.ip are dimensionless; therefore, the dimensional analysis of the 
space metric component gives us 



(9.5) 1 



2gM{f) 



M[g] 



L ' 

where [•] is the units of the variable. Rearranging this equation, we have 

m [ei = ^. 

In the same way that c = 1 makes the units of time equivalent to the units of length, 
setting Q = 1 equates units of length to units of mass. More specifically, let f^, be the radius 
measured in units of length, so the space coordinate in units of mass becomes 

(9.7) r- = |, 

where the space coordinate absorbs Newton's gravitational constant. This absorbtion is 
equivalent to setting ^ = 1 in (j9.ip . which is done in our simulations. For r in units of 
mass, we recover f^, in units of length by simply computing 



(9.8) 



r=K = Qr. 



To take this one step further, we also make t in units of mass related to in the units of 
length by 

(9.9) h = -t. 

c 

Therefore, by setting c = 1 and ^ = 1, we unify all the units (time, length, and mass) into 
one set of units. For our discussion, we choose the mass to be this unified set of units, where 
the units of length and time can be recovered with equations (|9.8p and (j9.9p , respectively. 

We also consider the units for the fluid variables {p,v). The velocity v is measured as 
a percentage to the speed of light and is thus unitless. With units of length equivalent to 
units of mass, the units for the density p in our simulation are 

M 1 

To recover the density in units of M/ , we compute 



(9.11) = ^p. 



To determine the value of Newton's gravitational constant, we perform dimensional 
analysis to recover the factors of c that have been suppressed by setting c = 1. Looking up 
the units, the Newton's true gravitational constant, denoted as Q, has units of 

7-3 

(9.12) r 



Mr2' 

Comparing (j9.6p to (j9.12p . the relationship between them must be 

(9.13) g = 

To obtain numerical values for Q and c, the units for the variables i^,, f^,, and M must 
be chosen. We choose to measure time in seconds (sec), the length scale in kilometers (km), 
and mass in solar masses Mq. This choice makes Newton's gravitational constant 



(9.14) 



G = 1.47664 km/M( 
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and the speed of light constant 

(9.15) c 3 X 10^ km/sec. 

With the abihty to recover the proper units, we interpret some of the numerical values 
of our simulated FRW/TOV model in Figure [7^ on a solar scale. To interpret the results 
on a solar scale, we choose the unit of mass to be represented by a solar mass (IM = IMq), 
where Mq = 1.989 x 10^° g is the mass of our sun. Using this scale, the range of values for 
the variables are 



(9.16) 0.09 Mq < M < 1.5 Mq, 

(9.17) 3 Mq < r<7 Mq, 

(9.18) 5.46 Mq < t<6.46M0, 

(9.19) 3.48 X 10"'' 1/M^ < P< 1-92 x 10"^ 1/M|,, 

(9.20) < V < 0.53. 

Using the constants ()9.14p and (j9.15p and equations (j9.8p and ()9.9p . we convert the mass 
to the familiar units of 

(9.21) 0.09 Mq < M < 1.5 Mq, 

(9.22) 4.43 km < u < 10.37 km, 

(9.23) 2.69 X 10~^ sec < U < 3.18 x 10"^ sec, 

(9.24) 1.08 X lO^'^Afo/km^ < p* < 5.95 x lO~''M0/km^ 

(9.25) km/sec < v < 1.59 x 10^ km/sec. 



Hence, our simulation is dealing with mass equivalent to 1.4 times our sun across a distance 
of 5.94 kilometers for a time frame of 4.9 microseconds. To give more meaning to the 
density, the average density of the sun is 7.04 X lO-^^M0/km^ so our simulation is dealing 
with densities of at least 1.53 x 10^^ times greater than the average density of the sun. In 
this coordinate system, matter is moving at speeds up to 1.59 x 10^ km/sec. 
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To consider galactic scales instead, one unit of mass is set equal to the mass of the 
Milky Way, 1.8 x 10^^ M©. The effect of this equivalency is a direct scale by the mass of the 
Milky Way of the solar scale results (|9.21|) - (|9.25p . except for the velocity. In this setting, 
the simulated values (|9.16p - (j9.20p become 



(9.26) 


1.62 X 10^° Mq 


< 


M < 2.7 X 10^^ Mq, 


(9.27) 


0.084 light-years 


< 


< 0.2 light-years, 


(9.28) 


56 days 


< 


t* < 66 days. 


(9.29) 


1.94 X lO'^Mo/km^ 


< 


< 1.07 X lO^Mo/km^ 


(9.30) 


km/sec 


< 


V < 1.59 X 10^ km/sec. 



Therefore, on the galactic scale, our simulation is dealing with mass equivalent to 1.4 times 
our Milky Way across a distance of 0.12 light-years in a time frame of 10 days. 
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APPENDIX A 

Programming Code for the Simulation 

This appendix provides an overview of all the source code used to perform the simula- 
tions throughout this thesis. A copy of this source code is provided on the attached CD, 
which is approximately 8,000 lines contained in 22 files. This code is written in C++, with a 
Microsoft C++ 5.0 compiler, using the OpenGL graphics library on the Windows operating 
system. Our goal here is to give the reader a sense of the organization of the code to enable 
the reader to find specific areas of interest within the code. 

To start, the main file (main.cpp) initializes the displayed window and the simulation 
itself, and contains the code specific to the Windows OS. This file also handles the interface 
between the operating system and the simulation. The rest of the code is segregated into 
different classes in accordance to an object-oriented programming paradigm. Each of these 
classes has two files associated with them, a header * . h file and an implementation file 
* . cpp. The header file contains the definition of the class and its associated variables and 
functions, and the implementation file contains the code for these functions. To reduce 
the amount of redundant code, we created a class hierarchy for the simulators, displayed 
in Figure [A) Classes on top of the hierarchy are base classes and the connected nodes 
represent children classes that inherit there functionality from the base class. We provide 
a brief explanation of the functionality of the classes in hierarchy: 

• The CSimulator class is a base class to provide a template in which all simulation 
classes are derived. It unifies the basic structure and interface of all of our simulator 
classes. 

• The CRsolver class is a base class to provide a template for a Riemann solver to the 
compressible Euler equations, both the relativistic and classical. The algorithms 
for finding the middle state and solving the Riemann problem, as discussed in 
Subsections 14.1.21 and 14.1.31 respectively, are contained in this class. Also, the 
code to display the Riemann problem simulator in Figure is in this class. 
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ClI CClasEulerSolw ]I> Cl] ^elEulerSolv ePl> 

Figure A.l. The object oriented structure of the classes in the simulation code 

• The CCosmoSim class contains the code to simulate and display the cosmology 
models of Chapters [6][8j This class contains the algorithm for the locally inertial 
Godunov method with dynamical time dilation, defined in Chapter [H Since the 
code became quite large, we split the implementation file into two, the numerics 
in c cosmos im. cpp and the graphics in ccosmosim_graphics . cpp. 

• The CRelEulerSolver class is the Riemann solver for the relative compressible Euler 
equations. This class is derived from the CRsolver class, and it contains the specific 
equations tied to the relative compressible Euler equations, scribed in Section 14.11 

• The CClasEulerSolver class is the Riemann solver for the classical compressible 
Euler equations [Smo83) . This class was used to test our Riemann solver on 
a more simple set of equations to find errors in the algorithms contained in the 
CRsovler class. This class has no affect on the simulations in this paper and is left 
over code from the testing stage of our simulations. 

To handle the rendering of our simulations, we use a number of stand alone classes 
to serve specific graphical needs. We present a list of these classes along with a brief 
explanation of their utility: 

• The Color class is used to consolidate the 3 color components into one vector. This 
class allows for easy manipulation and use of colors in OpenGL. 

• The Font class gives us a method to display all the text used in our simulations. 

• The Interval class is a concise method of keeping track of intervals on the real line. 



The ColorArray class enables us to keep track of an array of colors. 
The ColorMap class uses the Interval and ColorArray classes to build a mapping 
of an interval on the real line to a set of discrete colors. In particular, this class is 
used to display the density in an array of different colors in Figure 14.41 
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